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Abstract 


The  Three-Dimensional  Finite  Difference  (3DFD)  computer  code 
is  compared  using  Absorption  Boundary  Conditions  (ABS) 
versus  Radiation  Boundary  Conditions  (RBC).  This  comparison 
is  made  when  the  3DFD  code  is  used  to  study  the  interaction 
of  lightning  with  an  aircraft.  The  3DFD  computer  code  is  a 
modified  version  of  Rymes '  3DFD.  The  aircraft  modeled  for 
the  paper  is  an  F-16  'Fighting  Falcon'.  The  ABC  used 
simulates  an  infinite  free-space  by  setting  the  conductivity 
of  the  boundary  space  to  that  of  distilled  water,  to 
"absorb"  the  outgoing  electromagnetic  waves.  The  RBC 
simulates  free-space  by  assigning  the  boundary  fields  to  a 
previously  calculated  value.  The  value  is  calculated  with  a 
parabolic  interpolation  of  three  previous  field  values, 
which  are  offset  in  space.  Therefore,  the  calculated  value 
is  also  extrapolated  to  account  for  the  time  delay  and 
position  change.  The  results  of  incorporating  RBC  were 
dramatic.  The  ten  locations  sampled  for  the  test  showed 
marked  improvement  in  the  waveforms  when  using  RBC's. 
Depending  on  the  purpose  of  the  analysis,  this  improved 
waveform  output  may  be  overshadowed  by  the  25%  increase  in 
CPU  time  that  is  needed  for  the  more  sophisticated  RBC. 


Introduction 


1 . 


This  thesis  compares  Radiation  Boundary  Conditions 
(RBC)  with  Absorption  Boundary  Conditions  (ABC)  in  a  time- 
domain  three-dimensional  finite  difference  computer  code 
( 3DFD) .  In  this  application,  the  3DFD  computer  code  is  used 
to  analyze  lightning’s  electromagnetic  interaction  with  an 
F-16  aircraft.  The  major  tasks  accomplished  during  this 
thesis  effort  include; 


a.  The  F-16  'Fighting  Falcon*  aircraft  was 
electromagnetically  modeled,  and  implemented 
in  a  modified  version  of  the  Rymes ’  3DFD  code 
(16)  . 

b.  The  3DFD  code  with  its  original 
absorption  boundary  conditions  was  updated, 
corrected  and  run  for  a  typical  nose-to-tail 
lightning  strike  (Appendix  1)  (48).  The 
fields  calculated  were  recorded  at  10  monitor 
locations  on  the  F-16  aircraft  (Appendix  3). 

c.  The  3DFD  code  was  then  modified  to 
implement  the  radiation  boundary  conditions 
(Appendix  2)  (59). 

d.  The  3DFD  code  with  the  new  radiation 
boundary  conditions  was  run  for  the  same 
nose-to-tail  lightning  strike  on  the  F-16 


1 


geometric  model.  The  fields  calculated  were 
sampled  and  recorded  at  the  same  10  locations 
during  both  runs  (Appendix  3). 
e.  The  results  of  the  two  computer  runs  for 
the  first  two  microseconds  were  compared 
(Chapter  6).  The  comparisons  were  based  on 
one  complete  run  with  the  absorption  boundary 
conditions  and  one  run  with  the  radiation 
boundary  conditions.  The  comparisons  were 
made  for  the  same  10  sample  points  on  the  F- 
16  (Appendix  3). 

The  F-16  was  geometrically  modeled  and  implemented  into 
computer  code  in  a  manner  suitable  for  analysis  by  a 
modified  version  of  the  Rymes '  3DFD  computer  code  (48;  16). 
The  Rymes'  3DFD  code  was  modified  and  used  by  ILt  Hebert  at 
the  Air  Force  V/right  Aeronautical  Laboratories  Atmospheric 
Electricity  Hazards  Group  ( AFWAL/FIESL) ,  Wright-Patterson 
AFB,  Ohio  (16).  Hebert  implemented  the  code  with  absorption 
boundary  conditions,  and  modeled  a  CV-580  aircraft  for  his 
study  (16).  During  this  thesis  effort  radiation  boundary 
conditions  similar  to  those  discussed  in  a  paper  by  Kunz  and 
Lee  were  developed,  encoded,  and  written  into  the  modified 
3DFD  code  (27).  The  computer  codes  were  run  for  a  direct 
lightning  strike,  nose-to-tail ,  on  an  F-16  aircraft.  The 
results  are  compared  and  presented  in  this  thesis. 
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Lightning 

Lightning  brings  to  mind  a  flash  of  light  in  the  sky,  a 
loud  clapping  sound  and  rumbling  thunder.  This  flash  of 
light  and  the  associated  sounds  are  the  discharging  and 
neutralization  of  the  atmosphere's  large  charge  centers  from 
one  cloud  to  another,  or  from  clouds  to  the  earth  (61:1). 
Lightning  in  this  text  is  to  be  thought  of  as  a  high  current 
electric  discharge  which  has  a  path  length  measured  in 
kilometers  (61:1).  Aircraft  in  the  presence  of  these  highly 
charged  areas  can  become  part  of  the  high  current  channel 
(35:90-1).  In  fact,  Mazur  believes  the  aircraft  itself 
triggers  lightning  (35:90-2).  Shaeffer  agrees  with  this, 
except  with  the  qualification  that  a  highly  charged  air 
mass,  one  which  is  conducive  to  a  lightning  discharge,  must 
already  exist  before  the  aircraft  can  trigger  the  lightning 
discharge  (50:67). 


Lightning  and  Aircraft 

Lightning  poses  a  possible  catastrophic  threat  to  any 
aircraft  (ll:i).  Lightning  strikes  account  for  more  than 
one-half  of  the  total  USAF  weather-related  aircraft  mishaps 
(36:vi).  Since  little  can  be  done  to  prevent  an  aircraft 
from  being  struck  by  lightning,  aircraft  protection  for  the 
critical  components  is  of  utmost  importance  (45).  Critical 


components,  such  as  flight  control  computers,  fuel  systems, 
and  structurally  critical  airframe/propulsion  units,  are 
vulnerable  to  electrical  disruption  and/or  mechanical 
failure  (45).  As  electronic  complexity  has  increased,  the 
size  and  weight  have  decreased.  However,  the  probability  of 
lightning  strikes  has  remained  relatively  constant  (40:8). 

It  is  interesting  to  note  that  Pierce  concludes,  using 
basically  the  same  logic,  that  the  lightning  hazard  to 
aircraft  operation  is  increasing  (42:17).  Therefore,  the 
importance  of  characterizing  lightning  and  its  effects  is  of 
growing  importance,  as  indicated  by  the  number  of  articles 
written  on  the  subject  (8;  23;  31;  35;  36;  40;  42;  43;  46; 
47;  49;  50).  An  aircraft,  such  as  the  F-16,  with  a  "fly-by¬ 
wire"  flight  control  system,  were  it  unprotected,  could  be 
rendered  inoperable  by  a  lightning  strike  (1:1). 

The  starting  point  in  characterizing  lightning's 
effects  with  a  computer  program  is  to  model  the  aircraft's 
outer  skin.  Studying  the  propagation  of  the  electromagnetic 
energy  from  the  entry  to  the  exit  point,  as  it  redistributes 
on  the  surface  of  the  aircraft's  fuselage,  is  the  first  goal 
(34:2,16).  Determining  this  redistribution  can  lead  to 
coupling  the  surface  currents  to  the  interior  components. 

The  final  goal  is  understanding  lightning's  effect  on  the 
sensitive  electronic  components,  and  then  protecting 
sensitive  components  from  these  effects.  An  accurate 
computer  code  could  be  a  valuable  tool  in  designing 


protection  against  lightning's  effects.  It  would  allow 
protection  to  be  designed  and  tested  before  the  production 
stage  of  new  aircraft  development. 

All-metal  aircraft  are  shielded  from  many  of 
lightning's  effects,  as  the  all-metal  skin  presents  a 
Faraday  cage  type  enclosure  for  the  electronics  bay  (19;  45; 
52;  58).  Recent  efforts  to  use  the  new,  lighter,  advanced 
composite  and  thermoplastic  materials  in  aircraft  reduces 
this  shielding  protection  (10:1).  As  early  as  1952,  Burkley 
was  looking  at  plywood,  plastics,  laminated  fiberglass,  and 
other  non-conductive  materials  to  be  used  as  aircraft 
structural  components  (5).  Burkley  was  studying  the  effects 
of,  and  the  necessary  protection  needed  from,  the  high 
cu  rent  surges  associated  with  lightning  strikes  on  these 
non-conductive  materials  (5).  Examples  and  techniques  of 
various  protection  schemes  can  be  found  in  numerous  sources. 
Sommer  of  Boeing  Company,  Weinstock  of  McDonnell  Aircraft 
Company,  and  Fisher  of  General  Electric  are  just  three  of 
the  many  authors  publishing  in  this  area  of  aircraft 
lightning  protection  (51;  63;  14). 


Electromagnetic  computer  codes  allow  the  engineer  to 
characterize  and  study  the  electromagnetic  interaction  with 
aircraft.  The  engineer  can  then  design  protection  against 
those  effects  which  cause  electronic  upset  or  damage.  As 
previously  stated,  the  use  of  advanced  composite  materials 
and  sensitive  microelectronics  in  more  modern  aircraft  makes 
them  more  susceptible  to  the  effects  of  lightning  strikes. 
This  increased  susceptability  demands  more  accurate 
electromagnetic  modeling  (46:220).  This  modeling  is  to  be 
used  in  the  development  of  lightning  protection.  Accurate 
models,  when  used  in  design  for  lightning  protection,  allow 
the  use  of  optimal  lightning  protection.  This  protection  is 
optimal  in  that  excessive  lightning  protection  would  add 
weight  that  is  not  in  proportion  to  the  risk  reduction 
gained  (63:34-1). 

While  accurate,  these  models  must  also  be  flexible 
enough  to  handle  new  lightning  channel  models.  A  recent 
inflight  lightning  characterization  program  has  shown  that 
intracloud  lightning  attachment  is  even  more  severe  than 
previously  estimated,  with  faster  risetimes  and  higher 
charge  transfers  (19:1).  This  demonstrates  the  need  for 
flexible  analysis  systems  which  can  be  easily  modified  to 
incorporate  new  findings  without  completely  reworking  each 


Many  electromagnetic  computer  codes  are  available  to 
approximate  a  solution  to  Maxwell's  equations  (2;  13;  28; 


34;  60).  One  type  of  code  solves  Maxwell's  equations  in  the 
integral  form  based  upon  Harrington's  Method  of  Moments, 
such  as  Auckland's  study  of  an  F-14  (15;  4).  Another  type 
of  code  solves  Maxwell's  equations  in  their  differential 
form  using  the  Finite  Difference  Method  of  approximation, 
such  as  Holland's  "THREDE"  (20;  21;  22).  The  time-domain 
finite  difference  electromagnetic  code  approximation  of 
Maxwell's  equations  was  used  during  this  thesis  effort.  In 
a  study  by  Longmire  of  Mission  Research  Corporation,  it  is 
stated  that  the  finite  difference  method  is  the  fastest  way 
of  solving  Maxwell's  equations  (33s2341).  And,  in  a  recent 
article,  Mei  states  that  he  believes  the  finite  difference 
method  is  the  most  adaptable  to  the  virtual  memory  systems 
of  minicomputers  which  are  in  growing  use  in  industry 
(37:1145) . 


Background 

In  the  study  of  lightning's  interaction  with  aircraft, 
it  is  often  necessary  to  have  a  computer  model  for 
analytical  comparisons.  Hebert  has  successfully  modified  a 
finite  difference  electromagnetic  code  written  by  M.D.  Rymes 
of  Electromagnetic  Applications,  Inc.,  to  analyze  a  CV-580 
lightning  strike  aircraft  test-bed  (48;  16).  Other 


modifications  to  Rymes  3DFD  code  and  an  expanded  user's 
guide  are  detailed  in  Hebert's  report  (16).  This  modified 
3DFD  code  presently  uses  absorption  boundary  conditions 
(Appendix  1),  which  cause  the  boundaries  to  artificially 
reflect  unabsorbed  electromagnetic  energy.  This  gives  less 
than  desired  accuracy,  as  the  reflected  waves  return  to  the 
aircraft's  surface  at  later  time  steps  (55:626). 

3DFD  is  a  finite  difference  formulation  of  the  time- 
domain  electromagnetic-field  problem.  This  code's  major 
advantage  is  that  it  does  not  require  the  memory  or  the  time 
needed  to  invert  the  matrices  encountered  in  the  MOM  codes 
(33:2340).  As  Mur  states  it,  though,  the  finite  difference 
method  has  as  its  major  problem  the  limited  size  of  the 
problem  space  (41:377).  It  must  be  remembered  that  it  is 
the  problem  space  size  that  is  the  problem,  not  the  size  of 
the  object  inside.  For  example,  an  aircraft  as  large  as  a 
B-52  has  been  studied  using  3DFD  (21).  A  Cartesian  finite 
difference  method  uses  a  rectangular  problem  space,  totally 
enclosing  the  object  to  be  studied  (e.g.,  aircraft).  The 
object  to  be  studied  is  defined  in  the  problem  space  by 
assigning  values  of  permittivity  and  conductivity  to  each 
component-  of  the  total  electric  field  (E-field),  which 
describes  the  geometry  of  the  object.  This  geometry  has  the 
restriction  of  being  described  by  a  limited  number  of 
predetermined  rectangular  shapes,  which  make  up  the  problem 
space  in  the  Cartesian  coordinate  system. 


The  computer  code  solves  Maxwell's  time  dependent  curl 
equations,  which  in  turn  solve  the  boundary  conditions  for 
the  object  in  a  "natural  way"  (62:397).  The  code  progresses 
through  the  problem  space  in  time  steps  using  an  algorithm 
originally  developed  by  Yee  (64).  Since  an  infinite  problem 
space  cannot  be  defined  on  the  computer,  difficulties  arise 
when  the  propagating  wave  reaches  the  problem  space 
boundary.  These  boundaries  cause  reflections  unless  they 
are  modified  to  account  for  the  ideal  analytical  situation 
of  free  space.  Thus,  additional  algorithms  are  needed  to 
account  for  the  radiation  conditions.  Taylor  was  the  first 
to  implement  Yee ' s  algorithm.  He  used  absorption  boundary 
conditions  to  account  for  the  previously-mentioned 
reflections  at  the  problem  space  boundaries  (59:585). 

Yee  originally  started  with  "hard"  lattice  truncation 
(another  way  of  expressing  the  boundary  condition)  (64). 

Hard  lattice  truncation  is  defined  as  forcing  the  outside 
boundary  of  the  problem  space  to  be  a  perfect  conductor. 

This  is  done  by  assigning  the  tangential  E-field  the  value 
of  zero  at  the  boundary  (55:626).  This  is  also  known  as 
"tin  can"  boundary  conditions,  as  the  problem  space  is 
totally  enclosed  by  perfectly  conducting  metal  boundaries 
(16:22) . 

The  two  methods  that  are  examined  in  this  thesis  are 
both  referred  to  as  either  "soft"  lattice  truncation 
methods,  or  "soft"  boundary  conditions.  These  are  the 


absorption  boundary  condition  and  the  radiation  boundary 
condition.  The  absorption  boundary  condition  is 
accomplished  by  assigning  increasing  values  of  conductivity 
to  the  cells  of  the  problem  space  near  the  boundaries, 
thereby  effectively  chipping  away  at  the  E-field,  a  small 
amount  at  a  time  (55:626).  The  other  method,  the  radiation 
boundary  condition,  reduces  the  magnitude  of  the  E-field  by 
a  factor  of  l/r,  which  is  the  radiation  condition 
characteristic  of  an  electromagnetic  wave  propagating  in 
free-space  (39:41).  Here  'r'  is  the  radial  distance  from 
the  origin  of  a  centrally  located  coordinate  system. 

The  finite  difference  codes  have  been  found  to  agree, 
within  1  dB  and  1  lattice  cell,  with  known  analytical  and 
experimental  quantities  (54:202).  Due  to  the  accuracy  and 
efficiency  of  Yee's  basic  algorithm,  Taflove,  Kunz, 
Umashankar,  Merewether,  Fisher,  Mei,  and  others  have 
published  many  papers  on  combined  methods,  hybrids,  curved 
surfaces,  and  other  modifications,  making  this  a  very 
competitive  method  for  a  numerical  solution  of  Maxwell's 
equations  (14;  26;  37;  38;  56;  62). 


Problem 


First,  an  electromagnetic  model  of  an  F-16  aircraft  had 
to  be  implemented  on  the  computer  for  a  new  study  using  the 
present  modified  code.  During  the  computer  runs  using  the 
absorption  boundary  conditions  (ABC),  sample  data  from 
several  points  was  stored  for  further  study.  This  required 
the  correct  geometrical  modeling  of  the  aircraft  into  the 
three-dimensional  finite  difference  problem  space.  Next, 
the  code  was  modified  to  incorporate  the  radiation  boundary 
conditions  (which  were  believed  to  limit  previous 
reflections).  Then  the  code  was  run  with  the  same  F-16 
model,  sampling  the  same  points  as  before.  The  final  task 
was  to  compare  the  results  of  both  computer  runs  and  analyze 
the  findings. 


Scope 


This  study  was  limited  to  developing  a  geometrically 
correct  electromagnetic  model  of  the  F-16  using  the 
subroutine  AIRPLN  from  an  AFWAL/FIESL  technical  report  by 
ILt  Hebert  (16).  The  subroutine  AIRPLN,  modeling  an  F-16, 
was  run  with  the  present  modified  computer  code  for  a  nose- 
to-tail  lightning  strike.  Samples  were  recorded  at  10 
locations  (7  H-fields  and  3  E-fields)  for  predetermined 
orientations.  The  code  was  then  changed  to  incorporate  the 


radiation  boundary  conditions  (RBC).  The  modified  version 
of  the  code  with  radiation  boundary  conditions  was  run  with 
the  same  subroutine  AIRPLN  .  The  same  10  sample  points  were 
recorded  and  the  results  were  analyzed. 

Assumptions 

The  assumptions  were  that  the  AFWAL/FIESL  technical 
report,  and  previous  studies  which  determined  that  RBC  would 
be  more  accurate  than  ABC,  were  correct  (16).  An  isotropic, 
linear,  and  homogeneous  medium  was  considered.  The  region 
of  interest  was  considered  source-free  except  for  the 
injected  lightning  strike. 


2 .  Theory 


The  three-dimensional  finite  difference  code  uses  the 
time-domain  differential  form  of  Maxwell's  equations 
(64:302).  The  E-field  and  H-field  magnitudes  in  each  of  the 
three  vector  directions  that  define  the  Cartesian  coordinate 
system  are  calculated  using  these  equations  (64).  The 
calculations  are  made  while  stepping  through  the  grid 
problem  space  (NxNxN)  at  a  particular  time.  It  is  necessary 
to  increment  the  time  step  after  each  complete  pass  through 
the  grid  space  calculating  the  E-field  or  H-field.  All 
units  are  in  the  MKS  system  unless  otherwise  specified.  The 
form  of  Maxwell's  equations  for  isotropic  material  used  here 
is  (24:361) 


3H 

Vi~  =  -VxE  (2.1) 

3t 

3E 

G —  =  V  x  H  -  aE  (2.2) 

3t 

V  .  (eE)  =  p  (2.3) 


V  .  (pH)  =  0  (2.4) 


where 


H  =  magnetic  field  (amps/meter) 

E  =  electric  field  (volts/meter) 
t  =  time  (seconds) 

U  =  (mu)  permeability  (henry/meter) 
e  =  (epsilon)  permittivity  (farad/meter)  _ 

p  =  (rho)  electric  charge  density  (coulomb/meter"^) 
a  =  (sigma)  conductivity  (mho/meter) 

V  X  =  curl  operator 

V  •  =  divergence  operator 

These  are  point  relationships  (J  =  oE).  In  the  finite 
difference  code,  the  term  'decentralizing  mesh'  is  often 
used  to  refer  to  the  gridding  of  the  problem  space.  This  is 
not  to  be  confused  with  the  mesh  relationship  of  the 
integral  form  of  Maxwell's  equations. 

A  necessary  process  in  understanding  the  finite 
difference  code  is  to  develop  the  three-dimensional  finite 
difference  form  of  Maxwell's  equations.  This  development 
will  be  using  a  central  differencing  estimation  for  a 
derivative  (32:163).  The  first  operation  in  the  process  is 
to  put  Maxwell's  equations  in  a  finite  difference  form. 

Then  one  must  both  understand  where  the  fields  are  located 
in  the  problem  space,  and  comprehend  the  makeup  of  the 
problem  space  itself.  In  particular,  the  single  point 
reference  system  references  fields  located  on  three  sides  of 
the  rectangular  box.  A  point  in  the  three-dimensional  grid, 
and  the  associated  fields  are  shown  in  Figure  1.  Finally, 
one  must  put  all  of  this  together  and  express  the  three- 
dimensional  finite  difference  equations  in  a  form  useful  to 
algorithm  development.  A  source-free  region  is  considered 
for  the  development  of  the  algorithm. 


Figure  1.  Three-dimensional  Field  Representation 

Since  the  computer  memory  available  is  a  limited 
quantity,  some  boundary  must  be  placed  on  the  problem  space 
in  order  to  arrive  at  a  meaningful  solution.  Each  author 
has  his  own  style  of  implementing  boundary  conditions.  The 
two  compared  in  this  thesis  are  often  referred  to  as  "soft" 
boundary  conditions,  or  "soft"  lattice  truncation 
conditions.  The  first  to  be  used  is  the  absorption  boundary 
condition  which  Rymes  implemented  in  3DFD  (48).  The  second 
boundary  condition  is  the  radiation  boundary  condition  used 
by  Kunz  ( 29) . 


Finite  Difference  Form  of  Maxwell  *  s  Equations 

The  finite  difference  form  of  Maxwell's  equations  uses 
Eq  (2.1)  thru  Eq  (2.4),  manipulated  in  the  Cartesian 
coordinate  system,  to  arrive  at  a  convenient  algorithmic 
form  for  programming. 

In  the  Cartesian  coordinate  system,  the  electric  and 
magnetic  field  vectors  are  as  found  in  Thiele's  notation 
(53:11) ; 


A 


+  A  a  +  A  a 
y  y  z  z 


(2.5 


where  'A'  is  equal  to  the  magnitude  of  the  E-field  or 
H-field  component  and  'a'  is  the  unit  vector  in  the  x,  y,  o 
z  direction. 


The  curl  of  A  in  the  Cartesian  coordinate  system  is 
expressed  as  follows  (24:178); 


curl  A  =  V  X  A  =  det 


(2.6 


When  the  determinant  (det)  of  Eq  (2.6)  is  found,  it  is 


expressed  as  (59:556): 


(3 A  3A\^  /a a  3A  /3A  3A  \ 

— - +( — - +( — ^ - (2.7) 

3y  3z/  *  \  3z  3x/  ^  \  3x  3y/ 


In  general,  'A'  may  also  be  a  function  of  time. 


Substituting  E  or  H  for  A  into  Eq  (2.7),  replacing  the 
curls  in  Eq  (2.1)  and  Eq  (2.2)  with  the  results  from  Eq 
(2.7),  and  lastly  separating  the  vector  components  yields: 


3E^  3E, 


(2.8a) 


p  = 

3t 


+  !!5 

3z  dx 


(2.8b) 


3E  3E 

-  -J:  + 

3  X  3y 


(2.8c) 


+  o  E  = 
X 


3H  3H 

_ Z 


(2.9a) 


+  0  E  = 


3  Z  3  X 


(2.9b) 


+  oE  = 
z 


3H  3H 


3  X  3y 


(2.9c) 
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These  equations  are  now  in  a  form  convenient  for  use  in  the 
algorithm  of  the  finite  difference  computer  code. 


The  derivatives  in  Eq  (2.8)  and  Eq  (2.9)  are  replaced 
by  a  finite  difference  approximation.  The  exact  definition 
of  a  derivative  using  forward  differencing  is  (25:289; 

7:297)  : 

df(x)  f(x+h)  -  f(x)  ~  f(x+h)  -  f(x) 

f'(x)  =  -  =  lim  -  =  -  (2.10) 

dx  h— >0  h  h 

The  finite  differencing  approximation  used  in  this  work  is 
the  central  differencing  approximation  defined  by  (7:298): 


df(x)  ~  f(x+h/2)  -  f(x-h/2) 
dx  h 


and  is  illustrated  graphically  in  Figure  2. 


(2.11) 


(a)  Central 


(b)  Forward 


Figure  2.  Differencing  Points 

By  using  a  Taylor  series  approximation  of  Eq  (2.11)  and 
Eq  (2.12),  the  errors  can  be  compared  (32;  16).  The  result 
is  that  the  central  differencing  approximation  is  a  second 
order  approximation  as  opposed  to  forward  differencing, 
which  is  a  first  order  approximation  (32:297,  298).  Both  of 
the  codes  (Appendix  A  and  Appendix  B)  were  run  using  the 
central  differencing  approximations  for  each  derivative. 


Field  Locations 


The  points  (i,j,k)  which  describe  the  location  of  the 
fields  in  each  cell  space  are  reference  points  on  a 
decentralizing  mesh  scheme.  This  decentralizing  mesh  scheme 
is  best  described  by  first  considering  the  one-dimensional 
decentralizing  mesh  (Figure  3).  Note  that  the  E-fields  and 
H-fields  are  not  co-located.  This  displacement  keeps  one 
from  knowing  the  E-field  and  H— field  components  at  the  same 
point  in  space.  The  fields  are  only  known  at  1/2  the 
differential  distance  ( /ix/2,  Ay/2,  Az/2)  in  any  particular 
coordinate  system  direction.  Also  note  that  the  E-field  and 
H-field  are  not  known  at  the  same  time,  but  that  they  are 
separated  by  1/2  the  time  step  (At/2). 


The  Three-Dimensional  Finite  Difference  Equations  can 
be  derived  by  combining  the  differential  form  of  Maxwell's 
equations,  the  definition  of  curl  in  the  Cartesian 
coordinate  system,  and  the  central  differencing 
approximation  to  a  derivative. 

Magnetic-Field  Algorithm  Development .  An  example  in 
terms  of  would  combine  Eq  (2.8a)  and  Eq  (2.11)  at  a  point 


specified  as: 


and  yield: 


{ (i-l/2)  Ax,  j  Ay,kAz,  (n-1/2)  At) 


(i-1/2) AX, jAy»kAz,nAt]  -  { (i-l/2) Ax, jAy ,kAz, (n-1 ) At } 


E  {(i-l/2) AX, j Ay, (k+l/2)Az, (n-1/2) At) 


E{ (i-l/2) AX, j Ay , (k-1/2) Az, (n-1/2) At ) 


E^I (i-l/2) AX, ( j +1/2) Ay ,kAz, (n-1/2) At ) 


(2.12) 


E{(i-l/2)Ax,(j-l/2)Ay,kAz,(n-l/2)At) 

2 


In  using  this  equation,  an  approximation  for  can  be 
obtained  by  knowing  E  and  E  at  1/2  a  space  increment  on 


either  side  of  at  At/2  earlier  and  one  At  earlier. 

Figure  4  illustrates  these  relationships.  Taflove  and 
Umashankar  have  shown  that  this  approximation  is  within  the 
accuracy  of  computer  calculations  (54;  62). 


y 


Figure  4.  Location  of  Difference  Fields 

Note;  The  origin  is  centered  at 

{ (i-1/2) Ax, jAy.kAz} ,  at  the  time  nAt. 


Next,  the  E-field  can  be  determined  in  a  similar  manner 
At/2  later  at  1/2  a  space  increment  in  the  x  direction.  The 
program  steps  through  the  decentralizing  mesh  in  half-space 
increments  and  half-time  increments  until  the  entire  grid 
space  is  covered  (NxNxN). 
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At 


I  ^ 


E^{(i+l/2)Ax, ( j-1/2) Ay.kAz, (n-l/2)At) 


E^{(i-l/2)Ax, ( j-l/2)Ay,kAz, (n-l/2)At} 


E^{iAx,(j-l/2)Ay, (k+1/2) Az, (n-l/2) At ) 


E^{iAx, ( j -1/2) Ay, (k-l/2)Az, (n-l/2) At) 


at  point  location  {iAx, ( j-1/2 ) Ay ,kAz, (n-l/2 ) At }  and 


(2.13) 


H  {i Ax, j Ay , (k-l/2 ) Az, nAt }  -  H  {i Ax, jAy , (k-1/ 2) Az , (n-1 ) At ) 
z  z 


at  point  location  {iAx, jAy , (k-l/2) az, (n-l/2) At ] . 


Electric  Field  Algorithm  Development .  Operating  on  the 
E~field  in  Eq  (2.9a)  at  grid  point 

{ (i-1 ) AX, ( j-l/2) Ay , {k-l/2) Az,n At ) 

E^  is  first  needed  at  t  =  (n+l/2)At.  Using  the  average  of 
E^  at  At/2,  before  and  after  t. 


E  {(n+l/2)At)  +  E  {(n-l/2)At} 

E^(nAt)  =  — - - -  (2.15) 

2 


substituting  in  the  space  coordinates,  one  finds 

E^{ (i-1) AX, (j-l/2) Ay, (k-l/2) Az,nAt ) 

E^{(i-1)AX, (j-l/2) Ay, (k-l/2)Az, (n+1/2) At] 

2 

(2.16) 

E  {(i-1) Ax, (j-l/2) Ay, (k-l/2) Az, (n-l/2) At] 

+  _± - 

2 


and  substituting  into  Eq  (2.9a),  then  expanding  the  right 
side  as  was  done  in  Eq  (2.8),  the  result  is 


(e/At  +  o/2) ( { { i-1) Ax, { j-1/2) Ay# (k-1/2) AZ, (n+l/2) At ) ) 

-  (E/At  -  a/2)(E^{(i-l)Ax,{j-l/2)Ay,(k-l/2)Az,{n-l/2)At}) 


{ ( i-1 ) Ax , j  Ay , (k-l/2 ) Az, nAt } 

Ay 

{ ( i-1 ) Ax, ( j-1 ) Ay , (k-l/2) Az, nAt } 
Ay 

Hy{(i-1) Ax, (j-l/2) Ay,kAz,nAt} 

Az 

Hy{ (i-1) AX, ( j-l/2) Ay, (k-1) Az,nAtl 


(2.17) 


Now  sigma,  the  conductivity,  is  in  the  equation.  This 
IS  the  area,  that  if  one  were  dealing  with  something  other 
than  a  perfect  conductor,  changes  in  the  conductivity  could 
be  made.  In  the  data  base,  one  can  express  the  conductivity 
at  each  grid  point.  In  the  program  it  is  averaged  in  this 
manner : 


Ox  =  0  l(i-l) AX, (j-l/2) Ay, (k-l/2) Az) 

o{(i-l/2) AX, (j-l/2) Ay, (k-l/2) Az) 
2 


(2.18) 


0{(i-3/2)  AX,  (  j-3/2)  Ay,  (k-l/2)  Az) 


And  epsilon,  the  permittivity  is 


=  e{ (i-l)Ax, ( j-l/2)AY, (k-l/2)Az}  (2.19) 

Similar  operations  for  at 

{ (i-l/2)Ax, ( j-l)Ay, (k-l/2)AZ,nAtl 

yields 

(e/At  +  o/2)(Ey{(i-l/2)Ax,(j-l)Ay,(k-l/2)Az,(n+l/2)At)) 

-  (e/At  -  o/2)(Ey{(i-l/2)Ax,(j-l)Ay,(k-l/2)Az,(n-l/2)At}) 

(i-1/2) Ax, ( j-1) Ay,kAz,nAt} 

Az 

Hx{ (i-1/2) AX, ( j-1) Ay, (k-1) Az,nAt) 

AZ 

(2.20) 

H^{iAX,  (j-1)  Ay,  (k-1/ 2)  Az,nAt) 

AX 

H  {(i-1)  AX,  (  j-1)  Ay,  (k-1/ 2)  AZ,nAt) 

+  5 - 

AX 

and  for  E  at  { ( i-l/ 2 ) ax , ( j-1/ 2 ) Ay , (k-1 ) Az , n At ) 


and 


E^{i, j,k, (n+l/2)At}  =  E^{i, j,k, (n-l/2) At} 


^e/At  -  0/2'' 
^e/At  +  0/2) 


,(e/At  +  o/2)y 


( j-*-l/2)  ,k,nAt} 
<  Ay 


H  {i, ( j-1/2) ,k,nAt} 


Ay 

Hy(i, j, (k+1/2) ,nAt) 
Az 


(2.23) 


H  {i, j, (k-l/2),nAt} 
+  _i - - ^ - 


Az 


In  this  manner,  the  can  be  calculated  knowing  the  H-field 
At/2  earlier  at  the  four  adjacent  points  in  the  orthogonal 
plane,  Ay/ 2  and  Az/2  away,  and  the  E^  At/2  earlier  at  the 
same  point.  A  similar  statement  can  be  made  about  H^,  Eq 
(2.23)  . 


Thus,  at  every  time  step  (At/2),  a  complete  set  of 
either  E  or  H  fields  is  calculated  for  the  entire  problem 
space.  This  process  is  continued  in  time  until  a  steady 
state  response  is  reached. 

For  computational  stability.  At  must  satisfy  the 
Courant  stability  condition  (27:329;  64:303), 
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At  < 


(2.24) 


V ( Ax“^+Ay“^+Az“^) ^ 

where  is  the  maximum  velocity  of  propagation  in  the 

medium.  Eq.(2.24)  means  that  the  increment  of  the  time 
steps  must  be  smaller  than  the  time  it  takes  the  wave  to 
travel.  This  keeps  the  magnitude  of  the  field  finite 
differences  small,  so  that  the  solution  is  stable.  In  other 
words,  the  finite  differences  are  small  because  the  passage 
of  time  was  small  enough  to  allow  only  small  changes  in 
magnitudes . 


.  Boundary  Conditions 


Boundary  conditions  are  necessary  because  space  in  gen¬ 
eral  is  infinite,  but  the  computer  has  to  have  a  finite  num¬ 
ber  of  values  to  make  its  calculations.  An  infinite  space 
is  simulated  by  truncating  the  problem  space  with  carefully 
chosen  boundary  conditions.  Examining  Eq  (2.22)  and  (2.23), 
it  can  be  shown  (as  in  Hebert's  development  of  the  propaga¬ 
tion  in  one  direction)  that  some  type  of  boundary  condition 
is  absolutely  necessary  for  the  finite  difference  code  (16). 
Without  boundary  conditions,  the  fields  at  the  problem  space 
boundary  would  be  undefined  and  would  not  allow  the  calcula¬ 
tion  of  the  next  field  when  stepping  through  the  space.  The 
end  result  is  a  decreasing  data  base  of  past  fields  from 
which  to  calculate  the  subsequent  solutions. 

There  are  numerous  boundary  conditions  imposed  by 
different  authors  in  the  algorithms  of  the  finite  difference 
codes.  The  first  author  that  developed  the  finite 
difference  algorithm,  Yee,  simply  forced  the  outer  problem 
space  boundary  to  be  a  perfect  conductor  (64).  This  causes 
'noise*  in  the  solution,  due  to  the  reflections  off  the 
boundary  (55:626).  Others  use  a  combination  of  Yee's, 

"hard"  boundary  conditions  and  "soft"  ones,  which  reduce  the 
reflection  problem  (55:626).  The  two  "soft"  boundary 
conditions  considered  are  the  absorption  and  the  radiation 
boundary  conditions. 


Absorption  Boundary  Condition 


and  time  for  calculations  (55:626).  This  particular 
boundary  condition  was  used  in  Rymes '  code  (48). 


Radiation  Boundary  Condition 

The  radiation  boundary  condition  uses  a  far-field 
approximation  and  the  radiation  condition  to  emulate  an 
ideal  infinite  problem  space.  The  general  form  of  the 
radiation  condition  was  originally  introduced  by  Merewether 
(38:41).  Bayless  and  Turkel  have  shown  that  the  radiation 
boundary  condition  is  a  very  valid  mathematical  techniaue 
(4).  The  form  of  the  radiation  boundary  condition  is:  some 
function  of  time  f(t-r/c)  divided  by  radial  distance,  r, 
from  the  center  of  the  problem  space  (38:41). 


E 


f(t-r/c) 


r 


(3.1) 


where 


f  =  causal  vector  function 
t  =  time 

c  =  speed  of  light 

r  =  large  distance  from  the  center  of  the  test  object 


The  function,  'f',  used  by  Merewether  will  be  the  one 
incorporated  into  the  modified  Rymes'  code  (38:42).  The 
fields  at  the  boundary  are  found  by  parabolic  interpolation 


in  time,  of  the  fields  that  are  one  cell  in  from  the  outer 


boundary,  at  nAt,  (n-l)At  and  (n+l)At.  This  is  expressed  as 

E  at  the  outer  boundary  point  (i,j,k),  maximum  j -plane 
z 

boundary,  in  the  y-direction  as: 


where 
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(3.2) 


and  the  x,y,z  position  coordinates  for  the  outer  boundaries 
are  such  that  the  coordinate  system's  origin  is  in  the 
center  of  the  problem  space.  The  'c'  here  is  taken  to  be 
the  speed  of  light  in  free  space,  and  At  is  the  time  elapsed 
since  last  grid  position.  Similarly,  E^  can  be  expressed 
for  this  boundary.  Two  components,  for  each  of  the 
remaining  five  sides  of  the  problem  space,  are  also 
expressed  in  the  same  way.  Eq  (3.2)  works  for  a  constant 


cell  size  problem  space.  More  elaborate  expressions  are 
required  for  an  expanding  cell  size  grid  (20:2419). 

Implementing  Eq.(3.2)  into  code  turned  out  to  be  a 
major  task.  In  fact,  half  the  time  spent  on  this  thesis  was 
devoted  to  simply  understanding  this  radiation  boundary 
condition.  Researching  the  radiation  boundary  condition's 
origin,  and  then  implementing  it  into  the  3DFD  computer 
code,  was  a  much  larger  task  than  originally  estimated. 
References  were  not  found  which  specifically  covered  a  non- 
uniform  problem  space.  A  non-uniform  problem  space  in  this 
text  is  one  with  three  different  size  sides  to  each 


element's  rectangular  box.  A  major  clue  was  found  in  a 
report  by  Merewether  and  Fisher  (39:74).  Looking  at 
Eq.(3.2),  the  expression  for  9  contains  a  '1'  and  a  cAt. 

The  expression  directly  coded  into  FORTRAN  gave  disastrous 
results . 

This  direct  programing  caused  the  calculated  field 
magnitudes  to  blow  up.  In  looking  into  the  problem,  it 
became  obvious  that  cAt  had  to  be  different  for  each 
direction  considered.  This  was  taken  care  of  by  realizing 
that  cAt  was  the  distance  traveled  in  the  increment  of  time 
in  a  particular  direction.  This  could  simply  be  replaced  by 
the  space  increment  dimension  over  which  the  derivative  was 
being  considered.  Note  that  the  '1*  now  needs  to  represent 
an  integer  one  higher  than  the  quantity  being  subtracted 
away  in  the  expression  for  0  (7:120).  If  this  last 
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precaution  is  not  taken,  the  sign  of  0  varies  and  the 
parabolic  interpolation  is  incorrect.  Detailed  expressions 
can  be  found  in  Appendix  B. 

A  smoother  waveform  is  the  expected  result  when  using 
the  radiation  boundary  conditions  instead  of  the  absorption 
boundary  conditions.  The  reflected  waves  which  are  caused 
by  the  discontinuities  at  the  changes  in  problem  space 
conductivity,  are  expected  to  be  reduced.  Because  of  this, 
the  interaction  of  the  reflected  waves  should  be  greatly 
reduced  or  eliminated.  In  general,  the  reflective  waves  can 
be  thought  of  as  destructive  and  constructive  interference 
of  the  output/resulting  waveform.  This  interference  is 
believed  to  be  the  cause  of  the  abnormal  magnitude 
variations.  Therefore,  reducing  this  destructive  and 
constructive  interference  caused  by  the  reflected  waves 
should  result  in  a  smoother-appearing  output  waveform. 


L 


The  F-16  'Fighting  Falcon'  was  chosen  as  the  aircraft 
to  be  modeled  in  this  thesis.  The  F-16  was  chosen  because 
of  the  current  lightning  susceptibility  testing  being 
conducted.  The  thesis  sponsor,  AFWAL/FIESL,  was  performing 
lightning  strike  tests  on  the  new  LANTIRN  navigation  pod 
installed  on  a  F-16  (9).  Another  reason  for  the  choice  of 
the  F-16  was  due  to  its  being  the  first  operational  fighter 
with  a  "fly-by-wire"  flight  control  system.  The  high 
interest  in  the  survivability  of  the  "fly-by-wire"  system 
after  a  lightning  strike  caused  many  studies  to  be 
performed,  supplemented  with  many  on-going  r>f forts  to 
protect  this  system  (6).  A  large  data  base  was  therefore 
available  for  follow-on  comparisons  and  studies.  First,  the 
problem  space  will  be  described  in  general.  Then  the  choices 
which  were  made  to  model  the  F-16  in  the  time-domain  finite 
difference  problem  space,  will  be  explained.  The  last  area 
to  be  covered  will  be  the  location  of  the  ten  sample  points. 


Problem  Space 


The  problem  space  is  made  up  of  a  rectangular  area 
subdivided  into  19,683  rectangular  cells.  The  positions 
within  this  problem  space  are  defined  in  terms  of  the 
Cartesian  coordinate  system.  This  particular  choice  of 
coordinate  systems  is  well-suited  to  the  finite  difference 
formulation  previously  mentioned.  The  space  is  divided  into 
rectangular  grids  referenced  by  integers  such  as 
( i , j ,k ) =( 1 , 1 , 1 ) .  The  particular  length  of  each  grid  is 
chosen  to  enhance  the  description  of  the  geometry  of  the 
subject  object.  Each  of  the  three  directions  is  divided 
into  the  same  number  of  grids,  although  this  is  not 
required.  The  best  way  to  describe  this  is  by  an  example. 
The  program  used  for  this  thesis  incorporates  a  27X27X27 
(=19,683)  grid  space.  In  meters,  the  x-direction  increment 
is  0.69  meters,  y-direction  is  0.22  meters,  and  the  z- 
direction  is  0.45  meters  (Figure  6).  Since  the  grids  remain 
a  uniform  size  in  a  given  direction,  this  makes  the  overall 
physical  size  of  the  grid  space  18.63,  5.94,  and  12.15 
meters  respectively,  in  the  x,y,z  directions. 


z 


Figure  6.  Basic  Grid  Block 

The  constraint  of  memory  size  and  processing  time  are 
the  main  limiting  factors  for  the  number  of  grid  divisions 
used  (e.g.,  27X27X27).  Each  of  the  19,683  grid  points  has 
six  fields  associated  with  it.  After  one  complete  run  of 
the  program  through  the  problem  space,  118,098  fields  have 
been  calculated.  This  run  through  the  problem  space  occurs 
at  every  time  increment  (At/2).  The  time  step  in  this 
thesis  was  0.366667  nanoseconds.  A  valid  run  may  include  up 
to  two  microseconds  of  data.  More  time  results  in  large 
numerical  errors  (13:50).  Two  microseconds  corresponds  to 
5455  complete  cycles  through  the  problem  space.  Clearly, 
this  is  already  amounting  to  large  amounts  of  memory  and 
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central  processor  (CPU)  time.  When  the  amount  of 
calculations  per  cycle  is  considered,  the  problem  becomes 
evident . 

The  finest  detail  needed  to  be  represented  on  the  model 
determines  the  increment  size  in  each  direction.  Complete 
descriptions  can  be  found  on  the  problem  space  in  most 
finite  difference  users'  manuals  (18;  20;  29;  48;  60). 

These  same  users '  manuals  will  contain  information  on 
implementing  the  most  efficient  dimensions  for  the  problem 
space  and  element  blocks.  The  27  cubed  problem  space  was 
used  in  comparing  both  the  radiation  and  absorption  boundary 
conditions.  The  aircraft  was  allowed  to  take  up  a  large 
portion  of  the  problem  space. 

The  remainder  of  the  problem  space  left  a  minimal 
number  of  rectangular  blocks  beyond  the  aircraft  dimensions. 
The  blocks  outside  the  aircraft's  dimensions  are  used  for 
the  boundary  conditions.  Using  a  minimum  number  of  blocks 
is  not  a  good  practice  in  general,  but  was  intentionally 
done  to  exercise  the  boundary  conditions  under  a  worst  case 
scenario.  In  particular,  only  three  blocks  were  used  to 
implement  the  boundary  conditions.  The  aircraft's  maximum 
dimensions  were  extended  into  an  aircraft  grid  space 
21X21X21,  centered  in  the  27X27X27  problem  space.  This 
leaves  the  6  blocks  (3  on  each  face  of  the  box)  in  each 
direction  for  use  in  implementing  the  boundary  conditions. 


Field  Locations 


Each  corner  of  a  rectangular  box  in  the  problem  space 
defines  a  set  of  fields  (Figure  7).  The  6  fields  described 
are  the  3  electric  fields  and  the  3  magnetic  fields  in  each 
vector  direction.  The  fields  are  only  defined  by  these 
points;  they  are  not  actually  located  there.  The  E-fields 
defined  by  point  (i,j,k)  are  perpendicular  and  centered  on 
the  three  sides  of  the  rectangular  box,  not  touching  the 
defining  point.  The  H-fields  lie  on  the  edges  of  these  same 
sides  most  distant  from  the  defining  point.  The  directional 
nature  of  these  6  fields  is  as  established  by  the  coordinate 
system  used  (Figure  7). 


Figure  7. 


X 

Field  Locations  with  the  Decentralizing  Grid 
System 
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Cell  Size  Choices  Made  for  the  F-16 


The  F-16  is  a  very  smooth,  aerodynamically  designed 
aircraft.  This  aerodynamic  design  consists  of  many  curved 
surfaces  (Figure  8).  These  curved  surfaces  presented  a 
challenge  to  model  with  the  rectangular  blocks  of  the 
Cartesian  coordinate  system  in  the  problem  space.  The  F-16 
has  an  overall  length  of  15.09  meters,  width  of  9.45  meters, 
and  a  height  of  4.6  meters  (Figure  8)  (9;FO-3).  The  x 
coordinate  is  associated  with  the  length  of  the  aircraft 
(Figure  9).  The  x  direction  increment  of  0.69  meters  was 
chosen  as  a  compromise  to  better  model  the  geometrical 
structure  of  the  wings  in  the  xz  plane  (Figure  9) .  The  y 
coordinate  is  associated  with  the  height  of  the  aircraft 
(Figure  10).  The  y  direction  increment  of  0.22  meters  was 
selected  primarily  to  model  the  wing  root  area  of  the  F-16 
(Figure  10).  The  z  coordinate  is  associated  with  the  width 
of  the  aircraft  (Figure  11).  The  z  direction  increment  of 
0.45  meters  was  chosen  to  model  the  detail  of  the  speed 
brake  and  the  base  of  the  vertical  stabilizer.  The  result  of 


these  choices  was  a  good  geometrical  model  which  matches  the 
F-16's  predominant  details  (Figure  12). 
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Figure  11.  Gridding  of  the  F-16 


Figure  12.  Three-dimensional  Drawing  of  F-16  as  Blocks 

The  aircraft  is  described  to  the  computer  in  subroutine 
'AIRPLN'  (Appendix  A).  The  description  of  the  geometry  of 
the  aircraft  is  accomplished  by  defining  the  tangential  E- 


fields  of  the  surfaces.  The  building  of  the  model  (Figure 
13)  was  an  extremely  important  step  in  defining  the  location 
of  the  tangential  surface  fields.  Each  field  must  be 
assigned  by  grid  space  location  and  magnitude.  The 
tangential  surface  E-fields,  surface  normal  H-fields  and  all 
interior  fields  are  set  equal  to  zero  in  the  case  of  a 
perfectly  conducting  surface.  Refer  to  Appendixes  A  and  B 
for  listings  of  all  of  the  computer  codes  used,  including 
the  subroutine  AIRPLN.  A  useful  device  for  detailing  the 
field  locations  was  a  clear  plexiglas  cube,  with  the  field 
locations  accordingly  marked  such  as  Figure  1.  Another 
technique  used  in  modeling  the  F-16  was  to  take  the  wood  and 
metal  model  and  position  the  xz-planes,  or  layers,  of  the 
mock-up  on  a  scaled  grid  plane.  This  enabled  actually 
seeing  where  the  fields  were  physically  located.  The  coding 
of  the  subroutine  AIRPLN  and  understanding  the  field 
locations  was  greatly  simplified  with  this  technique. 


5.  Program  Details 


The  modified  Rymes '  code  was  reviewed  in  great  detail. 
Some  changes  were  made  to  the  modified  code  before  running 
the  F-16  with  the  absorption  boundary  conditions.  The  code 
was  then  modified  to  include  the  radiation  boundary 
conditions.  The  sensor  locations  were  selected  and  encoded 
in  subroutines  'EADV'  and  'HADV*.  Lastly,  the  source  was 
updated  to  reflect  parameters  recently  measured  in  flight 
(8;  17;  47;  49). 


Basic  Code 


The  assumptions  made  about  the  modified  Rymes’  3DFD 
code  were  basically  sound.  In  reviewing  the  modified  code, 
only  a  few  minor  omissions  and  corrections  v»ere  made.  The 
results  after  complete  computer  runs  on  the  F-16  model  were 
as  expected.  The  changes  made  to  the  modified  code  with  the 
absorption  boundary  conditions  enhanced  the  results 
(Appendix  C).  A  completely  updated  and  corrected  copy  of 
the  code  is  in  Appendix  A. 


Boundary  Conditions 


After  running  the  code  with  the  absorption  boundary 
condition,  the  radiation  boundary  condition  was  implemented. 
The  specific  listing  can  be  found  in  Appendix  B,  subroutine 
'RADBC  and  'HADV.  A  significant  amount  of  time  was  spent 
in  transforming  the  rather  simple-appearing  parabolic 
interpolation  of  (Eq  3.2)  into  Fortran  77  code.  The  basic 
problem  was  caused  by  a  lack  of  details  in  the  literature. 
The  missing  details  were  definitions  of  time  increments  in 
relation  to  the  three  different  dimensions  of  the  basic 
problem  space  block  (Figure  6).  In  particular,  the 
calculation  of  '0'  is  very  sensitive  to  the  definition  of 
c  At  ( Eq  3.2). 


Sample  Point  Locations 

Ten  locations  were  selected  as  typical  for  field  sensor 
layouts  on  the  aircraft.  Similar  locations  were  used  in  a 
CV-580  flight  test  conducted  by  AFWAL  (16).  Any  of  the 
three  orientations  of  either  the  E-field  or  the  H-field 
could  have  been  selected  at  any  point  in  the  problem  space 
for  monitoring.  The  selections  made  are  detailed  in 
Table  1. 
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TABLE  1 .  Ten  Sampled  Points 


!i 

i: 


Sensor 

# 

Type 

H  or  E 

Location 

Description 

Direction 

Component 

Coordinates 

(i/  j»)«) 

1 

H-f ield 

right  wing 

X 

(18,10,6) 

2 

H-field 

nose 

z 

(4,10,13) 

3 

H-field 

engine  feathers 

z 

(24,11,13) 

4 

H-field 

bottom  fuselage 

z 

(10,4,14) 

5 

H-field 

bottom  fuselage 

z 

(19,4,14) 

6 

H-field 

vertical  stab 

X 

(22,18,14) 

7 

H-field 

left  wing 

X 

(18,10,23) 

8 

E-field 

right  wing 

y 

(17,11,10) 

9 

E-field 

top  fuselage 

y 

(15,13,14) 

10 

E-field 

left  wing 

y 

(17,11,18) 

► 


H 


The  sensor  locations  are  shown  in  Figure  14.  Appendix  A, 
subroutine  'EADV'  and  'HADV'  show  the  code  used  to  sample 
these  points  at  each  time  step.  This  sampling  will  be  found 
under  'sample  points  of  interest’  at  the  end  of  both 


subroutines . 


Source  Function 


A  couple  of  the  changes  made  to  the  code  were  in  the 
source  function.  The  "transparent”  source  function  utilized 
in  the  modified  Rymes '  code  caused  some  unnatural  growth  in 
the  response  (19).  The  so-called  "transparent"  source  added 
the  scattered  fields  to  the  source  fields,  resulting  in  an 
ever-increasing  source  function.  Modifications  made  in  the 
code  appearing  in  Appendix  A  corrected  this  growth  and  gave 
much  better  results.  The  risetime  (10%  to  90%  of  the 
maximum  amplitude)  was  changed  to  100  nanoseconds  (Figure 
15).  This  was  found,  by  AFWAL  measurements,  to  be  more 
characteristic  of  cloud-to-cloud  lightning  strikes  (8:128). 
The  normalized  source  function  is  plotted  in  Figure  16  out 
to  the  50%  falloff  point,  which  occurs  at  20.4  microseconds. 

Numerous  programs,  or  researchers,  use  this  double 
exponential  as  the  source  function,  which  is  sometimes 
called  the  ' Bruce-Golde '  model  (8:36).  This  source 
function  emulates  a  lightning  strike’s  characteristic 
waveform  (19).  The  source  is  attached  to  the  aircraft  with 
a  surface  current  injection  technique,  as  described  by  Kunz 
(30:1423).  This  surface  current  injection  technique  is 
basically  accomplished  by  inserting  an  H-field  into  the  grid 
space  boundary  conditions.  The  insertion  of  the  H-field 
takes  place  around  the  point  that  the  lightning  attachment 
is  being  simulated.  The  H-field  is  approximated  from  the 
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relationship,  H  =  l/(2  r);  where  'I'  is  the  time-varying, 
analytical  current  model  for  the  lightning  channel  desired, 
and  ‘r*  is  the  radial  distance  from  the  point  of  desired 
injection  to  the  nearest  H-fields  surrounding  that  point 
(24:165) . 


AMPLITUDE 


6.  Analysis 


The  result  of  implementing  the  radiation  boundary 
conditions  appears  to  be  an  improvement  over  the  absorption 
boundary  conditions  results .  The  absorption  boundary 
condition  results  contai  i  many  fluctuations  in  amplitude 
that  can  be  attributed  to  the  boundary  conditions  (19;  55). 
When  the  waves  travel  and  hit  the  discontinuity  at  the 
element  rectangular  cells  which  have  increased  conductivity, 
they  reflect  back  with  a  discrete  magnitude.  These 
reflected  waves  cause  the  additional  amplitude  fluctuations 
not  found  with  the  radiation  boundary  condition.  All  ten 
sensors  showed  the  same  results;  a  general  smoothing  of  the 
amplitude  of  the  response  on  the  aircraft's  skin  (Appendix 
C) .  All  the  field  plots  (Appendix  C)  had  the  same  basic 
shape  before  and  after  the  boundary  conditions  were  changed. 
The  first  300  nanoseconds  of  each  fields  response  are 
plotted  in  this  chapter  (Figure  17  thru  Figure  26).  Each 
sensor's  response  is  plotted  on  the  same  graph  for  both 
absorption  boundary  conditions  and  radiation  boundary 
conditions.  The  absorption  boundary  condition  plots  are 
annotated  with  an  '*'.  The  smoother  response  of  the 
radiation  boundary  condition  resembles  the  MIE  solutions  as 
calculated  by  Holland  (20:206).  Note  that  the  disparity  in 
the  polarity  of  the  E-field  response  (sensor  8  &  10)  was 
possibly  caused  by  a  resonance  effect  in  the  wing  region. 
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Figure  20.  Sensor  Four,  H-field,  Forward  Fuselage  Bottom 


RMPS/METER 


Figure  21.  Sensor  Five,  H-field,  Rear  Fuselage  Bottom 
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Waveform  Appearance 


The  general  appearance  of  the  responses  from  the 
different  boundary  conditions  was  dramatic.  The  amplitude 
oscillation  frequency  of  the  radiation  boundary  condition 
response  was  less  distinct  than  that  of  the  absorption 
boundary  condition.  As  Figure  17  thru  Figure  26  display, 
the  radiation  condition  gave  a  much  smoother  response.  The 
overall  response  for  the  entire  two  microseconds  of  each  of 
the  computer  runs  was  consistent.  In  each  and  every  graph 
of  the  response,  the  radiation  boundary  condition  curve  was 
smoother  and  settled  to  a  steady  state  response  much  quicker 
than  that  of  the  absorption  boundary  condition. 

Note  also  that  the  amplitude,  in  general,  of  the 
radiation  boundary  curves  are  higher  than  those  of  the 
absorption  boundary  condition.  This  lower  absorption 
amplitude  was  due  in  part  to  the  source  itself  being 
absorbed.  That  is,  the  source  actually  goes  through  a 
couple  of  element  blocks  which  have  a  finite  conductivity. 
However,  a  small  portion  of  the  amplitude  differences  could 
be  due  to  constructive  and  destructive  interference  of  the 
reflected  and  surface  waves. 


The  absorption  boundary  condition's  H-field  response  of 
Figure  17  has  an  average  period  of  3.6  nanoseconds,  which 
corresponds  to  a  wavelength  of  1.08  meters.  The  radiation 
boundary  condition's  K-field  response  for  the  same  time  span 
is  different  (Figure  17).  It  has  a  average  period  of  7.2 
nanoseconds,  which  corresponds  to  a  wavelength  of  2.16 
meters.  Using  the  aircraft's  length  as  the  predominant 
resonance  structure,  the  wavelength  of  the  response 
calculated  for  radiation  boundary  conditions  was  a  quarter 
wavelength  multiple  of  the  aircraft's  length.  The 
absorption  boundary  condition  did  not  have  this  same 
relationship,  as  its  multiple  was  in  thirds  of  the 
aircraft's  length.  This  is  only  a  general  conclusion.  Time 
did  not  permit  taking  the  major  surface  dimensions  of  each 
sensor's  local  area,  and  correlating  the  calculated  field 
response's  frequency  to  the  major  local  surface  dimensions. 
But,  the  result  is  as  expected;  that  is,  the  major 
geometrical  object's  characteristic  wavelength  will 
dominate.  More  time  would  have  permitted  looking  at  the 
problem  space's  dominant  cavity  frequency,  and  its 
contribution  to  the  response  curve's  shape. 


Program  Considerations 

Implementing  the  radiation  boundary  conditions  has  many 
effects.  One  positive  effect  was  just  discussed,  that  of 
smoothing  the  response  waveform.  Other  effects,  such  as 
increased  storage  and  run  time,  are  not  desirable.  The 
amount  of  storage  for  running  the  program  of  Appendix  A  with 
its  absorption  boundary  condition  is  doubled  when  the 
radiation  boundary  condition  is  implemented.  An  additional 
363  lines  of  code  were  added  to  the  program,  to  implement 
the  radiation  boundary  condition.  The  worst  effect  was  the 
CPU  time  per  100  nanoseconds  of  program  time  increased  by 
25%  for  the  radiation  boundary  condition  over  the  absorption 
boundary  condition.  Table  2  contains  some  details  from  the 
CDC  Cyber  845  after  the  first  100  nanoseconds  of  data  was 
calculated.  No  claim  is  made  that  any  code  optimization  for 
run  time,  speed  or  memory  utilization  was  attempted. 

However  these  times  are  comparable  to  those  given  in 
Eriksen's  report  (13:50).  This  is  a  first  attempt  at  a 
comparison  on  the  same  code  with  only  the  boundary  condition 
changed.  The  similarity  of  the  plots  at  each  sensor  alone, 
support  the  validity  of  the  program  as  modified  with 
radiation  boundary  conditions. 


TABLE  2.  CDC  Cyber  845  Computer  Run  Data 


Absorption 

Boundary 

Condition 

Radiation 

Boundary 

Condition 

Percent 

Increase 

Time  elapsed 
within  program 

100  ns 

100  ns 

— 

SRU's  CDC  845 

8388.6 

11534.3 

37.5 

CPU  seconds 
execution  time 

3471.4 

4341.2 

25.1 

Number  cycles 
through  program 

272 

272 

- 

CPU  cost 
@  $320/hour 

$308.57 

$385.88 

25.1 

Cost  per  cycle 
(0.1833  ns) 

$1.13 

$1.42 

25.1 

Conclusion 


The  smoother  response  curves  of  the  radiation  boundary 
condition  are  quite  clear.  Assuming  this  smoothing  is  due 
to  ridding  the  solution  of  unwanted  resonances  by  changing 
the  boundary  conditions,  a  more  accurate  result  has  been 
accomplished.  But,  as  with  all  things,  more  study  needs  to 
be  done  before  real  confidence  can  be  placed  in  this 
computer  solution.  Checks  were  run  on  both  of  the  computer 
programs  to  verify  them.  Techniques,  such  as  running  the 
codes  without  sources,  or  without  an  object  in  the  problem 
space,  were  performed.  The  results  of  the  checks  of  the 


program  proved  that  it  worked  properly  in  those  aspects. 

But  more  experimentation  needs  to  be  performed.  On  the 
surface,  the  radiation  boundary  conditions  give  response 
curves  of  the  same  general  shape  as  found  in  lightning 
strike  characterization  papers.  This  comparison  is  made,  in 
general,  with  other  calculation  methods  of  response  curve 
generation.  A  lot  of  questions  yet  to  be  answered  have  to 
do  with  whether  or  not  the  increased  cost  of  this  program  is 
justified.  Can  the  analysis  be  accomplished,  with  the 
results  of  the  absorption  boundary  conditions,  to  the 
desired  accuracy?  Which  degree  of  accuracy  is  needed?  What 
is  the  decision  point  for  using  one  boundary  condition  over 
the  other? 

In  comparing  the  response  curves  side-by-side,  it  is 
noted  that  the  general  wave  patterns  for  each  sensor  are  the 
same.  Therefore,  implementing  the  radiation  boundary 
condition  did  not  alter  the  basic  wave  pattern.  TVnd  the 
radiation  boundary  condition  did  smooth  the  waveforms.  This 
tends  to  back  up  the  premise  that  the  spurious  reflections 
caused  by  the  absorption  boundary  condition's  changes  in 
conductivity,  are  seemingly  removed  by  the  radiation 
boundary  condition.  The  objective  of  ridding  the  output  of 
spurious  reflections  seems  to  have  been  accomplished. 


.  Summary  and  Recommendations 


The  major  objectives  of  this  study  have  been 
successfully  completed.  The  code  for  modeling  an  F-16  into 
the  Rymes ‘  3DFD  computer  code  was  written  and  implemented 
successfully.  The  3DFD  computer  code  was  run  with 
absorption  boundary  conditions  for  two  microseconds  of 
response  time.  The  3DFD  code  was  then  rewritten  to 
incorporate  the  radiation  boundary  conditions,  and  run  with 
the  radiation  boundary  conditions  for  two  microseconds  of 
response  time.  During  each  complete  run,  ten  points  on  the 
surface  of  the  simulated  aircraft  were  sampled.  These 
samples  are  in  the  form  of  magnetic  field  responses  and 
electric  field  responses. 

The  responses  for  each  sample  point  were  compared  for 
the  different  boundary  conditions.  The  radiation  boundary 
condition  produced  a  much  smoother  response,  apparently 
accomplishing  the  task  set  forth  by  the  sponsor, 
AFWAL/FIESL.  That  is,  the  radiation  boundary  conditions 
removed  many  of  the  unwanted  variations  in  the  amplitude. 
The  major  drawbacks  of  implementing  the  radiation  boundary 
conditions  over  the  absorption  boundary  conditions  were 
discussed.  The  25%  increase  in  the  cost  of  running  the 
radiation  boundary  condition  program  was  a  major 
disadvantage.  However,  the  overall  results  verified  that 


the  radiation  boundary  conditions  do  remove  unwanted  problem 
space  reflections.  The  sponsor  ( AFWAL/FIESL)  was  very 
pleased  with  the  results. 

Future  projects  to  parallel  this  effort  could  be 
numerous.  Comparisons  should  be  made  of  3DFD  output  data 

with  measured,  actual  inflight  lightning  strikes.  Another 
possibility,  would  be  exploring  in  greater  detail  the  result 
of  this  study  for  more  subtle  effects  of  the  different 
boundary  conditions  that  were  not  readily  apparent.  Or 
researchers  could  take  on  the  very  challenging  task  of 
optimizing  the  code  for  minimum  run  time  and  maximum 
efficiency  in  storage  use.  A  detailed  examination  should  be 
made  to  explain  why  the  E-fields  of  sensors  8  and  10,  seem 
to  have  opposite  polarities  for  the  different  boundary 
conditions.  An  interesting  technique  would  be  to  use  the 
aircraft  symmetry  in  such  a  manner  that  only  half  the 
problem  need  be  solved,  which,  in  theory,  would  reduce  the 
expense  by  half.  Researching  combinations  of  boundary 
conditions  would  also  be  of  great  interest.  One  possibility 
would  be  to  allow  a  small  amount  of  absorption  in  the 
outermost  cell,  in  combination  with  the  radiation  condition. 
Another  possibility  would  be  to  create  variations  in  the 
cell  dimensions  in  a  given  direction,  as  opposed  to  adding 
extra  grid  cells  (which  cause  cost  increases  or  loss  of 
details  in  the  modeling  of  the  object  being  studied). 


Many  studies  have  been  performed  in  the  area  of  Finite 
Difference  Codes,  as  attested  to  by  the  many  listings  in  the 
bibliography.  However,  few  real  comparisons  of  the  varied 
techniques  are  available  in  any  detail.  Much  can  be 
accomplished  with  this  code.  The  3DFD  is  very  flexible  and 
allows  cell-by-cell  assignment  of  material  properties.  This 
flexibility  makes  this  code  an  issue  of  possibly  great 
interest  in  the  study  of  composite  aircraft  materials  in 
many  electromagnetic  pulse  environments. 

Future  implementation  on  specialized  parallel 
processing  computers,  and/or  the  CRAY  super  computer,  would 
greatly  reduce  the  major  disadvantage  of  the  3DFD  code, 
which  is  the  CPU  runtime.  Of  course,  the  increased  cost  of 
the  use  of  these  newer  computers  must  be  weighed  against  the 
benefits.  The  last  suggestion  that  will  be  made  in  this 
thesis  is  that  an  interactive  modeling  system  be  developed 
for  electromagnetic  codes  in  general.  Modeling  an  aircraft 
is  a  very  time-consuming  task.  Cooperation  with  the 
aeronautical  engineers  and  their  detailed  drawing 
capabilities  would  open  new  frontiers  for  the 
electromagnetic  community.  The  ability  to  take  the 
aeronautical  engineer's  detailed  structural  drawing  and 
converting  them  into  electromagnetic  models  by  software, 
would  be  a  great  stride  forward  for  any  electromagnetic 
computer  code  utilization. 
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Appendix  A 


Updated  Modified  Rymes '  Code 


This  appendix  contains  the  Modified  Rymes'  Code,  with 
absorption  boundary  conditions,  written  by  Hebert  (16).  The 
code  listing  is  complete  as  implemented  and  runs  on  a  CDC 
Cyber  845.  Improvements  are  included  which  correct  minor 
errors  in  the  code  mechanics,  and  an  improved  source  is  also 
implemented.  The  code  compiles  on  the  Cyber  845  using  FTN5, 
and  runs  with  no  problem  when  suppressing  ANSI  errors.  The 
compile  time  is  approximately  2.603  CPU  seconds.  The  run 
time  is  3,472  CPU  seconds  for  100  nanoseconds  of  sample  time 
elapsed.  The  aircraft  modeled  in  this  code  is  an  F-16 
'Fighting  Falcon'.  This  code  will  run  as  listed  only  for 
the  first  150  nanoseconds.  For  continued  runs, 
modifications  must  be  made  as  annotated  in  the  program. 
Additionally  'TAPEl'  must  be  saved  from  the  previous  run. 


PROGRflM  ABC  ((]UTFlS,TflPEH]UTF15,TAPE2) 

C 

C  THE  HODIFIED  RYHES  3-D  CODE  FIE 
C 

COmON  /5ET1/  NNAX,EPST,SieO,C,TMJO,DELr 

COMMON  /SETE/  ICAN,ICM1,ICP1,JCAN,JCM1,JCP1,KCAN,KCM1,KCP1 

COMMON  /SET3/  DELX, DELY, OELZ 

COMMON  /TIME/  TEN,THN 

COMMON  /ARRAYS/  RHE7,27),B1(39),A2(27,27I,B2(39),A3I27,E7), 

4B3 (39) , A4 (27, 27) , B4 (39) , A5 (27, 27) ,  B5(39>, RE(27, 27) , BE (39) , 
W(27,27),B7(39),Aa(27,27),Ba(39),A9(27,27),B9(39), 
tAie (27, 27) , Bie(39) , A1 1 (27, 27) , B1 1 (39) 

DIMENSION  INDEX (190) 

CALL  0PENMS(2, INDEX, 190,0) 

C 

C  READ  INPUT  DATA 

C 

CALL  SETUP 
C 

C  ZERO  THE  FIELDS  INITIALLY. 

C  PROGRAM  USES  A  MASTER  STORAGE  FILE  TO  STORE  FiaDS  AND  CONDUCTIVITY 
C  OF  EACH  LATTICE  POINT.  INDEX  K  GOES  FROM  1  TO  27. 


c 

EX 

K 

c 

EY 

Ki27 

c 

EZ 

Ki54(27i2) 

c 

HX 

Kiai(27i3) 

c 

HY 

Kil0e(27i4) 

c 

HZ 

KM35(27i5) 

c 

SI6 

Kil62(27iG) 

C 

C(imi«f»iiffff  miff  iiiffiif  Ilf  fiiiifiimiiiifm 
CiiiiiiifiiiiDELETE  FROM  HERE  TO  STOP  DaETEiuiiiii 
C  TO  REMOVE  THE  ZEROING  OF  TAPEl 

C  FOR  CONTINUED  RUNS  BEYOND  THE 

C  INITIAL  TMAX  STORING  TIME 

DO  100  1=1,27 

DO  100  J=l,27 

100  R1(I,J)=0.0 

DO  200  K=l, 189 

200  CALL  URITMS(2,A1(1,1),768,K,0) 

C 

CiiifififiiiiSTOP  DaETE  60  THE  nW/TENiiiiiiiiiui 
Cifffffiifiififiiiiiiiiiiiifiiififiiiiiiiiffiififiif 

C  PRESET  TIC  CONDUCTIVITY  ARRAY 

C 

CALL  SIGSET 
C 

C  TIME  STEP  LOOP 

C 

DO  1000  N=1,NMAX 
CALL  APERT 
C 
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COMPUTE  TINE  CONSTANTS 


HHfMHitiFOR  CONTINUED  RUNS  BEYOND  TMflX«iM»f*»i( 
ADD  THE  STDPPIN6  TINE  FROM  PREVIOUS 
RUNS  TO  TUN  AND  TEN  PLUS  1/2  OF  DELTA  T 


EXAMPLE 

AFTER  FIRST  RUN 

THN=DaTi  (FLOAT  INI -I.  I)  ♦!.  501B3E-7 
TEN=DELTt  (FLOAT  (N) 5)  M.  50183E-7 

imiiiffiiKMtfiitiHiiiiitiiiiiiiiiiimifiiiiH 

THN=DaT»(aOAT(N)-1.0) 

TEN=DELTt (FLOAT (NI-. 51 

ADVANCE  nC  H-FIELDS  BY  Z-PLANES 

CALL  HADV 

ADVANCE  THE  E-FiaOS  BY  Z-PLANES 

CALL  EADV 
m  CONTINUE 

END  OF  TINE  LOOP 

CALL  CL0SNS(2I 

STOP 

END 

SUBROUTINE  SETUP 

CODE  USES  A  UNIFORM  GRID. 

COMMON  /SETI/  NMAX,EPST,S1GO,C,TMUO,D£LT 

COMMON  /5ET2/  ICAN,  ICNl,  ICP1,JCAN,JCMI,JCP1,KCAN,KCM1,KCPI 

COMMON  /SET3/DELX,DELY,DELZ 

INPUT  DATA  PROVIDED  BY  USER. 

DATA  C,  EPSO, EPSR,  XMUO,  5I(^/3. BED,  B.  B542E-I2, 1. 0, 1. 2SEGE-B,  I .  E-4/ 
DATA  ICPl,JCPl,KCPI/27,*:7.,27/ 

DATA  DELI,DaY,DELZ/.69,.22,.4S/ 

DATA  TMAX/1.5E-7/ 

C 

ICAN=1CP1-1 

1CM1=ICAN-1 

JCRN=JCP1-1 

JCM1»JCAN-1 

KCAN=KCPt-l 


KCMl=KCm-l 


VERIFY  THAT  ICPl, JCPI.KCPl  ARE  ODD. 

IF ( ICAN. ED. 2> ( ICAN/2) . AND.  JCAN.  EO.  2* ( JCAN/E) . AND. 
tKCAN.E0.2i(KCAN/2M60  TO  1 
PRINT  TIE  NUMBER  OF  GRID  POINTS  MUST  BE  ODD  * 

STOP 

CONTINUE 

BELT  IS  TIE  TIME  INCREMENT  (IT  SATISFIES  COURAMT  CONDITION). 

NMAK  IS  THE  TOTAL  NUMBER  OF  TIME  INCREMENTS  FiaDS  UOULD  BE  ADVANCED. 

DEIT=AMIN1  (DELX,  DELY,  DaZ)  /  (2.  <0 

NMAX=TMAX/DaT 

EPST=EPSOiEPSR/DaT 

TMUO=DaT/XMUO 

RETURN 

END 

SUBROUTINE  SI6SET 

SI8(I,J)  AT  K1=K41G2  IS  THE  CONDUCTIVITY  AT  THE  POINT  (X(I),Y(J),Z(K)) 
X(I)=(I-.5)»DaX,  Y(J)=(J-.5l»DaY,  Z(K)=(K-.5)iDaZ 

COMMON  /ARRAYS/  SI6(27,27),B1(39) 

COMMON  /SET!/  NMAX,EPCT,SIGO,C,TMUa,DELT 

COMMON  /SET2/  ICAN, ICN1,ICP1,  JCAN, JCMI,JCPI,KCAN,KCM1,KCP1 

PRESET  THE  CONDUCTIVITY  TO  lE-4  EVERYWERE  IN  SPACE. 

DO  10  I-I,ICAN 
DO  10  J=1,JCAN 
SIG(I,J)=SIGO 

NON  URITE  THE  SI6  ARRAY  TO  WSS  STORAGE 

DO  20  K=1,KCAN 
Kl=Ktl62 

CAa  NRITMS(2,SIG(1,I),7G8,K1,1) 

THEN  ESTABLISH  FRa  SPACE  CONDITIONS 

ICM2=ICRN-2 
JCM2=JCAN-2 
KCM2=KCAN-2 
DO  30  1=3,  ICM2 
DO  30  J=3,JCM2 
SIG(I,J)=0.0 

NON  OVERWRITE  THIS  PORTION  OF  THE  SIG  ARRAY 
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K1=KM62 

cm  URITNS(2,S!G(1,1),76B,K1,I) 

RETURN 

EhS) 

SUBROUTINE  EflOV 

THIS  SUBROUTINE  ADVANCES  THE  E-FIELDS  TUO  Z-PLANES  AT  A  TIME. 

RECORDS  OF  H-FIELDS  FDR  TUO  SUCCESIVE  Z-PLANES  MUST  BE  KEPT  IN 
ORDER  TO  TAKE  DERIVATIVES  IN  THE  Z-DIRECTION. 

RECORDS  OF  SIS  FOR  TUD  SUCCESIVE  Z-PLANES  MUST  BE  KEPT  IN  ORDER  TO  TAKE 
AVERAGES  IN  THE  Z-DIRECTION. 

COMMON  /SET!/WAX,EPST,SIGO,C,TMUO,DELT 

COmON  /SETR/  ICAN,ICM1,ICP1,JCAN,JCMI,JCPI,KCAN,KCMI,KCP1 

COMMON  /SET3/  DELX,DaY,DELZ 

COMMON  /TIME/  TEN,THN 

COMMON  /ARRAYS/  EX(27,27),BU3S),EYI27,27l,B2(39»,EZ(e7,27), 
+BI(3S),HXA(27,27),B4(39),HXB(27,27I,B5(39),KYA(27,27),B6(39), 
4HYD(27,27),B7I39),HZAI27,27),BS(39),HZB(27,27),B9I39), 
4SI6AI27,27),B10(39),SI6DI27,27),B11(39) 

ADVANCE  E-FIELDS  BY  Z-PLANES 
READ  IN  FIRST  H-PLANES 

CALL  READMS(2,HXA(I,1),7BS,B2) 

CALL  READMS(2,HYA(l,l),76a,ie9) 

DO  FOR  ALL  Z-PLAfCS 

MN8=fl 

MN9=9 

MNie=ie 

DO  IMB  KL0QP=1,KCAN,2 

K=KLOOP 

KPI=K4l 

K2=K427 

K3=K45* 

K4=K4B2 

K5=K+109 

K6=K4I35 

R7=K4l62 

READ  IN  SECOND  H-PUVCS  AND  OLD  E-FIELDS 


CALL  READMS(2,EXn,I),7Ga,K) 
CALL  READMS(2,EY(I,1),768,K2) 
CALL  RERDMS(2,EZ(I,1),768,K3I 
CALL  READMS(2,HXB(1,1),768,K4) 
CALL  READMS(2,HYDn,II,788,KS) 
CALL  RERDMS(2,HZA(t,I),7B8,K6) 
CALL  READMS(2,S16AU,I),768,K7) 


IWj 


DO  S09  >l,JCflN 
JPl=J+l 

DO  200  1=2, ICAN 

SIG(Xe,Y,Z)=.St(SIG(X-DELX,Y,Z)4SIG()(,Y,ZI) 
GIG(X,Y0,Z)=.5>(&I6(X,Y-I]ELY,Z)4&16(X,Y,ZI» 
GIG  ( X,  Y,  ZO)  =.  S>  (&I6(  X,  Y,  Z-DELZ)  +GI6  (X,  Y,  Z) ) 

6I62=.25i(GI6fl(I,J)4&IGA(I-l,J)) 

fl=EPGT-GIG2 

B=1./(EPST4GIG2) 

EX(I,J)=(niEXII,J)t(HZA(I,JPi)-HZfl(I,J))/DELY 
4(HYB(I,J)-HYn(l,J))/DELZ)iB 
CONTINUE 
DO  300  1=1, ICAN 
IP1=I41 

DO  300  J=2,JCAN 

GI62=.25t(GlGA(I,J-l)4SIGA(I,J)) 

A=EPST-G1G2 

B=l./(EPST4Gie2) 

EY(I,J)=(AiEYII,J)4(HXD(I,J)-HXA(I,J))/DElZ  • 
4<HZA(IP1,J)-NZA(I,J))/DELX»B 
CONTINUE 

URITE  NEW  EX,EY  PLANEG  TO  NASS  STORAGE 

CALL  URITI1S(2,EX(l,i),768,K,n 
CAa  WRITNS(2,EY(i,l),7£S,K2,l) 

IF  (K.E0.n6a  TO  450 
DO  400  1=1,  ICAN 
IPi=l4i 

DO  400  J=i,JCAN 
JP1=J41 

SI62=.25i(SIGA(I,J)4SIGD(I,J)) 

R=EPST-S1G2 

B=1./IEPST4SI62) 

EZ(I,J)=(AiEZ(I,J)  +  (HYA(IPl,J)-HYA(I,J))/DaX 
4lHXA(I,JPl)-HXA(I,JI)/DELYIiD 
CONTINUE 

URITE  NEW  EZ  PLANE  TO  NASS  STORAGE 

CALL  URITNS(2,EZ(l,l),768,K3,n 

SANPLE  POINTS  OF  INTEREST (ODD  K  PLANES) 

FORMAT  (I4,2(1X,)PIE11.4)) 

FORMAT  (14,2(IX,1PIE11.4),2(/)) 


CONTINUE 


CJ  O  CJ 


KP1=«+1 
Ke=K2tl 
K3=K3+1 
K4=K4+1 
K5=K5+I 
K6=»<6+1 
R7=K7+l 

READ  IN  SECOND  OLD  E-PLANES  AND  NEXT  H-PLANES  OFER  A  PLANE 

CALL  READHS(2,EX(l,i),768,K) 

CALL  READNS(a,EY(l,l),768,R2» 

CALL  READI1S(2,EZ(!,1),7E8,R3) 

CALL  READ)1S(2,HXA(l,n,7B8,RA) 

CALL  READnS(2,HYA(l,l),7SB,KS) 

CALL  READMSI2,H2e(l,l),7S8,K6) 

CALL  READMS(',SIGB(i,l),7S8,R7) 

00  £88  J=i,JCAN 
JP1=JM 

DO  £88  1=2, ICAN 

SI62=.25*(S16B{I,J)4SIGB(I-1,J)I 
A=EPST-SIG2 
B=1./(EPST4S1G2) 

EX(I,J)=(A»EX(I,J)4(HJB(l,JPi)-HZB(I,Jn/I)eLY  - 
♦(HYAII,J)-HYB(I,J))/DELZ)«B 
£88  CONTINUE 

DO  788  1=1, ICAN 
IP1=I41 

DO  788  J=2,JCAN 

SI62=.2St(SI6B(I,M)4SIGB(I,J)) 

A=EPST-SIG2 

B=I./(EPST4SI62) 

EY(I,J)=(AiEY(I,J)  +  (HXA(I,J)-HXB(I,J))/DaZ  - 
4(HZBIIP1,J)-HZB(I,J))/DELX)*B 
788  CONTINUE 

DO  888  1=1, ICAN 
IPI=l4l 

DO  888  J=1,JCAN 
JP1=J41 

SI62=.2St(SIGB(I,J)4SIGA(I,J)) 

A=€PST-SIG2 

B=1./1EPST4SIB2) 

EZ(I,J)  =  (A»EZ(I,J)4(HYB(IPl,J)-UYB(I,J)»/DaX  - 
4  (HXB n,  JPl )  -HXB 1 1,  J) )  /DELY)  iB 
888  CONTINUE 
C 

C  WRITE  NEN  EX,EY,EZ  PL^ES  TO  NASS  STORAGE 
C 

CALL  URITN5(2,EX(1,1),7£8,R,!) 

CALL  URITnS(2,EY(l,l),7£B,K2,n 
CALL  WRITNS(2,EZ(l,i),7£6,R3,l) 
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:  SRHPLE  roiKTS  OF  INTEREST  (EVEN  K  PLANES) 

<% 

IF  (K.EO.I0)  WRITE  (4,111)  >KB,TEN,EY(17,1I) 

IF  (K.EQ.IB)  WRITE  (4,112)  mi0,TEN,EY(17,ll) 

C 

IF  (K.EQ.14)  T)€N 

WRITE  (4,111)  m9,TEN,EY(iS,13) 

ELSE 
END  IF 
C 

1000  (XMTINUE 
C 

C  ZERO  TIC  aECTRIC  FiaO  WITHIN  THE  )CTH.  (»(PO)fNTS  OF  AIRCRAFT. 

C 

CALL  AIRPLN 

RETURN 

END 

SUBROUTINE  HADV 
C 

C  THIS  SUBROUTDt  ADVAtCES  THE  H-FiaOS  TWO  Z-PLANES  AT  A  TINE. 

C  RECORDS  OF  E-FiaOS  FOR  TWO  SUCCESIVE  PLANES  NUST  BE  KEPT  IN  ORDER  TO  TAKE 
C  DERIVATIVES  IN  THE  Z-DIRECTION. 

C 

COmON  /SETl/  NNAX,EPST,SieO,C,TNUO,DELT 

COItlON  /SET2/  ICAN,ICHl,lCPl,JCAN,JCNi,JCPl,KCAN,KCNl,KCPl 

CONMON  /SET3/  DELX,DaY,DaZ 

COWIQN  /TINE/  TEN,THN 

CONHON  /arrays/  E1(A(27,27),B1(3)),EXB(27,27),62(39), 

♦EYA(27, 27) , B3 (39) , EVB(27, 27) , B4 (39) ,EZA(27, 27) ,B5(39) , 

4EZB(27, 27), B6(39),HT  (27,27), B7(39),KY(27,27),B8(39), 
tHZ(27,27),B9(39) 

C 

C  DO  FOR  ALL  Z-PLANES 

»I(I=1 
NN2=2 
)f(3=3 
NN4=4 
NN5=5 
NN&=6 
NN7=7 
C 

DO  900  KL00P=1,KCAN,2 

K=KLOOP 

KNl=K-l 

K2=K*27 

K3=K»54 

K4=Ktfll 

K5=KM0B 

K6=Kfl35 

CALL  READNS(2,EXBI1,1),76B,K) 

CALL  READN5I2,EYB(I,1),768,K2) 

CALL  RERDH5(2,HZ(I,1),7E8,K6) 


IF  (K.EQ.1)GQ  TO  350 

CALL  READM5(»,EZB(t,l),7E8,K3) 

CALL  READNS(a,HX(l,l),76S,KM 
CALL  READMS(E,HY(1,1),7G8,KS) 

DO  200  J=2,JCAN 
DO  200  I»l, ICON 

HX(I,J)=HXLI,J)*THJOt((EYB(I,J)-EYO(I,Jn/DELZ  - 
4(EZB(I,J)-EZB(I,M))/DELV) 

200  CONTINUE 

DO  300  M,ICN1 
IPl=Itl 

DO  300  >1,JCAN 

HVLlPl,J)=HY(IPl,J)tTMUO»((EZB(IPl,J)-€ZBII,J))/DELX  - 
iLEXBdPl.JI-EXAdPi.JD/DaZ) 

300  CONTINUE 
C 

C  AOPROXINATES  H  FIELDS  IN  PLANE  PERPENDICULAR  TO  WIRE 
C  BY  H=I/2iPItR. 

C  SOURCE  FUNCTION, ODD  K  REFERENCE  PLANES,  HY 

C  CURRENT  PULSE  MOVING  IN  X  DIRECTION  ATTACIES  TO  TUE  NOSE  OF  THE  AIRPLANE 
C  AT  (X,Y,Z)=(B,ie,14),  EXITS  BY  THE  TAIL  AT  (X,Y,Z)«(22,1B,14). 

C  AT  TIME  T=0.0  PULSE  IS  LOCATED  AT  XO(S). 

C 

IFIK.NE.  15)60  TO  1 
DO  10  1=2,3 

TPULSE*THN-. 69» (FLOAT (I ) -3. I /C 
IF(TPULSE.LE.O.O)GO  TO  10 
C 

X*EXP(-3.5E04iTPULSE)-EXP(-2.625E7iTPl)LSE) 

C 

HY5=7.07E3»X 
HYd,9)=-HYS 
10  CONTINUE 

DO  30  1=25,26 

TPULSE=TH(-.  69f  (FLOAT  ( 1)  -3. 0)  /C 
IF(TPULSE.LE.O.OIGO  TO  30 
X=EXP (-3. 5E4iTPULSE) -EXP (-2. 625E7«TPULSEI 
HYS=7.07E3tX 
HYd,IO»=-HYS 
30  CONTINUE 

1  CONTINUE 

C 

C  WRITE  NEW  HX,HY  PLANES  TO  HASS  STORAGE 

C 

CALL  WRITMS(2,HX(I,1),768,K4,1) 

CALL  WRITMS(2,HYd,I),768,K5,l) 

350  CONTINUE 

DO  400  I=1,ICN1 
IP1=I*1 

DO  400  J=I,JCM1 
JPl=Jd 

HZ(IPl,JPl)=HZ(IPl,JPl)tTMUOt((EXB(IPl,JPl)-EXBdPI,J))/DELY  - 
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n  n  r> 


4(EYBIIPi,JPl)-EYB(I,JPI))/|]aX) 

400  CONTINUE 

URITE  NEU  HZ  PLANE  TO  MASS  STORAGE 

CALL  URITI1S(2,HZ(l,l),7Ba,KS,l) 

C 

C  SAMPLE  POINTS  OF  INTEREST  (ODD  K  PLANES) 

C 

IF  (K.EQ.I3)  WRITE  (4, III)  )te,Tm,HZ(4,)0) 

IF  (K.EQ.13)  URITE  (4, III)  MN3,T)K,HZ(24,1I) 

IF  (K.EQ.83)  URITE  (4,111)  )fl7,T)K,HX(lB,ie) 

C 

118  FORMAT  (I4,8(1X,1P1E11.4),8(/)) 

111  FORMAT  (I4,8(1X,1P1E11.4)) 

C 

K=Ktl 

KM1=K-1 

K8=K8+1 

K3=K3+1 

K4=K4+1 

K5=K5tl 

KS=KB4l 

C 

C  READ  IN  NEXT  PLANES  OF  E  AND  H 
C 

CALL  READn5(8,EXA(l,l),7Ea,KI 
CALI  READMS(8,EYRIl,l),7Ba,K8) 

CALL  READM5(8,EZA(l,ll,7Ga,K3) 

CALL  RF''.i:i«(8,HX(l,l),7Ga,K4) 

CAll  READM5(8,HY(l,l),76a,KS) 

CALL  READMS(8,HZ(1,1),7GB,K6) 

DO  S00  J=1,JCM1 
JP1=JM 

DO  800  1=1, ICAN 

HX(I,JPl)=HX(l,JPl)4TMUOi((EYAII,JPl)-EYB(I,JPl))/DELZ  - 
«(EZA(I,JP1)-EZR(I,J))/DELY) 

800  CONTINUE 

DO  700  I=1,ICM1 
IP1=I+1 

DO  700  M,JCAN 

HY(IPl,J)=«Y(lPl,J)*TMUOi((EZO(IPl,J)-EZA(I,J))/DELX  - 
«(EXA(IP1,J)-EXB(IPI,J))/DELZI 
700  CONTINUE 
C 

C  SOURCE  FUNCTION,  EVEN  K  REFERENCE  PLAIES,  HY 
C 

IF  (K.NE.  14)60  TO  2 
DO  80  1=8,3 

TPULSE=THM-.  89f  (FLOAT  ( I )  -3. )  /C 
IF  (TPULSE.LE.a.eiGO  TO  80 

X=EXP  ( -3. 5E4f  TPULSE)  -EXP  T-8.  B85E7»TPl)lSE) 


rj  o  r> 


HYS=7.07E3il( 

HY(I,g)»+HYS 
2»  CONTINUE 

DO  Al  1-25, % 

TPULSE=THM-.  69f  laOflT  ( 1)  -3. )  /C 
IF  (inJLSE.LE.e.eiGQ  TO  Ae 

X=EXPI-3. 5EAiTPULSE)-EXP(-2. 625E7fTPULSE» 

HYS=7.07E3»X 
HY(I,10)=+HYS 
A0  CONTINUE 

2  CONTINUE 

DO  800  I=I,ICN1 
IP1=I+I 

DO  000  M,JCN1 
JPl=Jtl 

HZ(IPl,JPl)=HZ(IPl,JPn*TNUOi(IEXn(IPl,JPn-EXA(IPl,J))/DELY  - 
♦  IEYfl(IPI,JPI)-EYA(I,JPn)/DaX» 

800  CONTINUE 

C 

C 

C  SOURCE  FUNCTION,  EVEN  K  REFERENCE  PLAf€,  HZ 

IF(K.NE.1A)G0  TO  A 
DO  80  1=2,3 

TPULSE=TIW-.  691  CROAT  ( 1)  -3. )  /C 
IF (TPULSE.LE. 0.0)00  TO  80 

X= (EXP  (-3.  SEAMPULSE)  -EXP  (-2. 625E7tTPaSE)  I 
HZS=I.AA7EAiX 
HZ(I,9)=-HZS 
HZ(I,10)=4HZS 
80  CONTINUE 

DO  00  1=25,28 

TPULSE=TH)t-.  83#  CROAT  ( I ) -3. ) /C 
IF(TPULSE.LE.O.O)GO  TO  00 
X=CEXPC-3.5EA»TPULSEI-EXPC-2.825E7*TPaSE)) 

HZS=I.AA7EA»X 
HZ(I,I0)=-MZS 
HZ(I,I1)=HZS 
00  CONTINUE 

A  CONTINUE 

CAR  WRITMSC2,HXCI,1),768,KA,I) 

CALL  HRIT«S(2,HYCI,l),760,K5,l) 

CALL  HRITHS(2,HZCI,l),780,K8,n 

SAMPLE  POINTS  OF  INTEREST  (EVEN  R  PLR)ES) 

IF  CK.E0.1A)  THEN 
MRITE  CA,1I1)  NNA,THN,HZC10,A) 

WRITE  CA,1II)  NN5,niC,HZCI9,A) 

WRITE  (A,l!l)  NN8,T)fC,HX(22,10) 

RSE 
END  IF 
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.  .'V" 


t  v  1 


it 

i. 


112  FiD(i,B)=e.e 
C  miN  FUSELAGE 

128  IF  (K.LT.12.0R.K.6T.1S)  GO  TO  126 

D0  121I»ie,19  I 

DO  121  J=S, 12 

121  Fij)(i,j)=e.e 

c 

C  FORUARD  FUSELAGE 

C  j 

DO  122  1=6,9 
DO  122  J=8, 10 

122  FLD(I,J)=8.I 
DO  123  1=8,9 
DO  123  J=1 1,12 

123  FLD(I,J)=8.0 
FLD(7, 111=8.8 

C  REAR  FUSELAGE 
DO  124  1=28,22 

124  aD(I,6)=e.0 
DO  12S  1=28,23 
DO  125  J=7, 12 

125  FLDn,J)=8.8  ^ 

126  IF  (K.NE.12.0R.K.NE.151  GO  TO  138 
DO  127  1=19,20 

D0  127J=3,4 

127  aD  n,J)=0.8 

FLD(28,5)=0.0 

138  IF  (K.LT.13.0R.K.6T.14)  GO  TO  135 
aO(24,9)=0.0 
aD(24,10)=8.8 

C  NOSE  AND  COCKPIT  : 

DO  131  1=4,5 

DO  131  J=8,9  , 

131  aD(I,J)=0.8  » 

aD(5, 181=0.8 

FLD(6, 111=8.0 
FLD (7, 121=0.0 
1)0  132  1=8,11 
DO  132  J=I3, 14 

132  aD(I,Jl=0.0  ‘ 

aD(12, 131=0.8 

FLD(13, 131=0.0 
aD(9, 151=0.0 
C 

C  RIGHT  VERTICAL  STABALIZER 

C 

135  IF  (K.NE. 131  GO  TO  140 

DO  136  1=19,23 

DO  136  3=13,14 

136  aDU,J)=0.0 
aD(17, 131=0.0 
aOllB.  131=0.0 


IF  (K.NE.U)  GO  TO  168 

DO  141  1:^18,23 

DO  U1  J=13,14 

FLD(I,J)=8.e 

FLO (17, 23) =8. 8 

DO  142  1^8,23 

DO  142  J=15,19 

FLD(I,J)«:8.8 

FLD(28,15)=8.8 

FlJ)l2e,  16)  =8. 8 

DO  143  1=22,23 

DO  143  3=28,22 

FLD(I,J)=8.8 

FLD(24,21)=8.8 

FLO (24, 22) =8. 8 

DO  144  1=23,24 

DO  144  J=23,24 

FLD(I,J)=8.8 

IF  (K.NE.16)  GO  TO  178 

DO  161  1=11,22 

DO  161  J=B,  18 

FLD(I,J)=8.8 

FLD(23,9)=8.8 

FlJ)(23,18)=e.8 

IF  (K.NE.17)  GO  TO  188 

WING 

DO  171  1=13,19 
DO  171  J=9,18 
aD(I,J)=8.8 

RIGHT  HORIZONTAL  STOBALIZER 

DO  172  1=21,24 
aD(I,9)=8.8 
IF  (K.NE.1B)  GO  TO  198 
DO  181  1=14,19 
DO  181  J=9, 18 
aD(I,J)=8.8 
DO  182  1=22,24 
aD(I,9)=e,8 

IF  (K.LT.19.0R.K.6T.28)  GO 

DO  191  1=15,19 

DO  191  J=9, 18 

aD(I,J)=8.8 

aO  (23, 9)  =0.8 

aO  (24, 9)  =0.8 

IF  (K.LT.21.0R.K.GT.22)  GO 

DO  211  1=16,19 

ao(i,9)=e.o 

IF  (K.NE.23)  GO  TO  248 

DO  231  1=17,19 

FLD(I,9)=8.8 


241  IF  (K.NE.24)  GO  TO  1000 
FLD(18,9)=0.0 
FLD(19,9)=0.0 

C  STORE  EX  FIELDS  BOCK  IN  TO  HASS  STORAGE 
1000  CALL  URITNS(2,FLD(1,1), 788, K,n 
C  ZERO  EY 

DO  2000  K=3,24 
L=K+27 

CALL  READNS  (2,FLD(I,1),768,U 
IF  IK.LT.7.0R.K.GT.I0)  GO  TO  1110 
DO  1081  1=14,19 
1081  FLD(I,1O)=0.O 

IF  (K.NE.9)  GO  TO  1100 
FLD(13,1O)=0.O 

1100  IF  (K.NE. 10)  60  TO  1110 
FLD(12,10)=0.O 
FLO(13,10)=0.O 

1110  IF  (K.NE. 11)  GO  TO  1120 

DO  nil  1=10,22 
DO  1111  J=9,10 
nil  aD(i,j)=0.o 

aD(23,lO)=0.O 
C  HAIN  FUSELAGE 

1120  IF  (K.LT.1'.0R.K.6T.1S)  GO  TO  113 

DO  1121  1=9,19 

DO  1121  J=G,12 

1121  aD(I,J)=0.0 

C  FUD  FUSELAGE 

DO  1122  1=5,8 
DO  1122  J=9,10 

1122  aD(I,J)=0.8 
DO  1123  1=7,8 
DO  1123  J=n,12 

1123  FLD(I,J)=0.O 
FLD(S,!1)=0.0 

C  REAR  FUSELAGE 

DO  1124  1=20,22 
DO  1124  J=7,12 

1124  FLD(I,J)=0.0 
DO  1125  J=0,12 

1125  FLDI23,J)=0.0 

1130  IF  (K.LT.13.0R.K.GT.14)  GO  TO  112 

FLO  124, 10) =0.0 

1126  IF  (K.NE.12.0R.K.)E.15)  GO  TO  113 

DO  1127  1=18,20 

DO  1127  J=4,5 

1127  FLD  n,J)=0.O 
FLD(20,6)=0.O 

1131  IF  (K.LT.13.0R.K.6T.14)  GO  TO  113 

FID(3,9)=0.0 

FLD(4,9)=0.0 

FLO(4,1O)=0.O 


FU)(S,ll)=8.e 
FLD(e,i8)=0.e 
DO  1132  h7,13 

1132  FlD(I,13)=e.e 
DO  1133  1=7,11 

1133  aD(I,U)=8.e 
RJ)(B,15)=«.0 
FLD(9, 151=0.0 

1135  IF  (K.NE. 131  GO  TO  1140 
DO  1136  1=16,23 

1136  FLD(1,13)=0.0 
DO  1137  1=10,23 

1137  FLD(I,14)=0.0 

1140  IF  (K.NE. 14)  GO  TO  1160 

DO  1141  1=22,23 

DO  1141  J=13,24 

1141  FLD(I,J)=0.0 
DO  1142  J=21,24 

1142  FLD(24,J)=0.0 
DO  1143  J=20,22 

1143  FLD(2t,J)=0.0 
DO  1144  J=17,19 
DO  1144  1=20,21 

1144  FlJ)(I,J)=0.0 
DO  1145  1=19,21 
DO  1145  J=1S,16 

1145  FLD(I,J)=0.0 
DO  1146  1=17,21 

1146  RDd, 141=0.0 
DO  1147  1=16,21 

1147  RDd,  131=0.0 

1160  IF  (K.NE. 16)  60  TO  1170 

DO  1161  1=10,22 

DO  1161  J=9,10 

1161  aD(I,Jl=0.0 
aD(23, 101=0.0 

C  HING 

1170  IF  (K.NE. 17)  GO  TO  1180 

DO  1171  1=12,19 

1171  FUld,  101=0.0 

1180  IF  (K.NE. 181  GO  TO  1190 

DO  1181  1=13,19 
1161  FLDd, 101=0.0 

1190  IF  (K.LT.19.OR.K.GT.20)  GO  TO  2000 
DO  1191  1=14,19 

1191  aDd, 101=0.0 

C  STORE  EY  BOCK  TO  MOSS  STORflGE 

2000  CALL  URlTNS(2,aD(l,  11,768,1,1) 

C  EZ  ZEROED 

DO  3000  K=3,24 
L=Kt54 

CALL  RE0DR5(2,FU)d,  11,766,1) 
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UIN6  TIP 

IF  (K.NE.4)  60  TO  SflSO 
DO  2041  1=17,19 

FLD{I,9)=0.I 
IF  (K.NE.S)  GO  TO  2060 
DO  2051  1*16, 19 
FIJ}(I,9)=0.0 
PLPNE 

IF  (K.LT.6.0R.K.6T.7)  6 
DO  2061  1=15,19 
FLD(I,9)=0.0 
UIN6 

IF  (K.NE.S)  GO  TO  2090 
DO  20B1  J=9, 10 
DO  2081  1=14,19 
FLD(1,J)=0.0 

RIGin  HAND  GTRBflLIZER 

DO  2082  1=22,24 
FLD(I,9)=0.I 
IF  (K.NE.9I  GO  TO  2100 
DO  2091  1=14,19 
DO  2091  J=9, 10 
FLDII,J)=0.0 
DO  2092  1=22,24 
aD(l,9)=0.0 
IF  (K.NE.10)  GO  TO  2110 
DO  2101  1=13,19 
DO  2101  J=9,18 
FU)II,J)=0.0 
DO  2102  1=21,24 
aD(I,9)=0.0 
IF  (K.NE.ll)  60  TO  2120 
DO  2111  1=12,24 
FLD(I,9)=0.I 
DO  2112  1=12,19 
FID(I,10)=0.0 

IF  (K.NE.12)  GO  TO  2130 
DO  2121  1=10,22 
DO  2121  J=e,10 
aDII,J)=0.0 
aD  (23, 9)  =0.0 
FID  (23, 10)  =0.0 

IF  (K.LT.13.0R.K.6T.I5) 
DO  2131  1=9,19 
DO  2131  J=5,12 
aD(I,J)=0.0 
DO  2132  1=5,8 
DO  2132  J=8,10 
aD(l,J)=0,0 
DO  2133  1=7,8 


DO  S133  J=ll,12 
FLD(l,J)=0.e 
FlDI6,ll)=e.8 
DO  2134  1=20,22 
DO  2134  J=6,12 
aD(I,J)=0.0 
DO  2135  J=7,12 
FLD  (23,31=0.0 
IF  (K.NE.14)  GO  TO  21G0 
FLO  (24, 9)  =0.0 
FID  (24, 10)  =0.0 
DO  2141  1=3,4 
DO  2141  3=8,9 
FLD(I,3)=0.0 
FLD(4,10)=0.0 
FID  (5, 111=0.0 
FLD (8, 12) =0.0 
DO  2142  1=7,11 
DO  2142  3=13,14 
aD(I,3)=0.0 
aD(12,13)=0.0 
FLD(!3,13)=0.0 
aD(S,lS)=0.0 
aD(9,15)=0.0 
DO  2143  1=16,23 
FLD(I,13)=0.0 
DO  2144  1=18,23 
aD(I,14)=0.0 
IF  (K.NE.16)  GO  TO  2170 
DO  2161  1=10,22 
DO  2161  3=8,10 
FLD(I,3)=0.0 
aO  (23, 9)  =0.0 
FU)(23, 10)  =0.0 

IF  (K.NE.17)  GO  TO  2180 
DO  2171  1=12,24 
FLO(I,9)=0.0 
DO  2172  1=12,19 
aD(l,10)=0.0 
IF  (K.)E.18)  GO  TO  2190 
DO  2181  1=13,19 
DO  2181  3=9,10 
aD(I,3)=0.0 
DO  2182  1=21,24 
aD(I,9)=0.0 

IF  (K.LT.19.OR.K.GT.20)  GO  TO  221 
DO  2191  1=14,19 
DO  2191  3=9,10 
aD(I,3)=0.0 
DO  2192  1=22,24 
aD(I,9)=0.0 

IF  (K.LT.21.0R.K.6T.22)  GO  TO  2230 


Appendix  B 


Rymes '  Code  with  Radiation  Boundary  Conditions 


This  appendix  contains  a  complete  listing  of  the 
modified  Rymes '  code  with  the  radiation  boundary  conditions 
implemented.  The  compile  time  is  3.722  CPU  seconds.  The 
run  time  for  the  first  100  nanoseconds  was  4342  CPU  seconds. 
The  F-16  is  modeled  in  subroutine  AIRPLN,  as  it  is  in  the 
previous  appendix.  The  major  changes  include  expansion  of 
the  memory  required  and  a  complete  rewrite  of  subroutine 
SIGSET,  now  renamed  RADBC.  The  code  listing  is  complete  as 
implemented,  and  runs  on  a  CDC  Cyber  845.  Improvements  are 
included  which  correct  minor  errors  in  code  mechanics,  and 
also  an  improved  source  is  implemented.  The  code  compiles 
on  the  Cyber  845  using  FTN5,  and  runs  with  no  problem  when 
suppressing  ANSI  errors.  Continued  runs  require  the  changes 
annotated  in  the  code  in  addition  to  saving  'TAPEl'  from  the 
previous  run. 


PROGRflN  RBC  (0UTF16,TnPE6=aJTFlC,TAPE  1) 

C 

C  THE  HODIFIED  RYIES  3-D  CODE  FIE 

CUIIfH 

C  TO  IMIUDE  RADIATION  BOUNDARY  CONDITIONS  DEVELOPED 
C  DEVELOPED  BY  NEREUETICR  AND  INCORPORATED  INTO  THIS 
C  CODE  BY  UILLIFORD  AND  (EBERT 

CtllHft 

c 

COmON  /SETl/  NWI, EAST, SIGO,C,TMUO, CELT 

COmON  /SET2/  ICAN,ICMl,ICPl,JCAN,JCMI,JCPl,KCAN,KCMi,KCPl 

COraiON  /SET3/  DELX,DELY,DaZ 

CON10N  /TINE/  TEN,THN 

COMMON  /ARRAYS/  A1(27,B7),B1(39(,A2(27,27),B2(39),A3(E7,B7), 

4B3 (39) , M (27, 27) , B4 (39) ,  AS (27, 27) , BS (39) ,  A6(27, 27) , B6 (39) , 

4A7 (27, 27) , B7 ( 39 ) , AS (27, 27) , 68 (39) , A9 (27, 27) , B9 (39) , 

4A10(27, 27) , B1B(39) ,  At  I (27, 27) , B1 1 (39) , A12(27,27) , B12 (39) , 

Ciiiiiti 

4A13(27,27),B13(39),A14(27,27),BU(39),A15(27,27),B15(39), 

4AlS(27,27),Blfi(39),A17(27,27),BI7(39),AIB(27,27),BlB(39), 

4A19(27,27),B19(39),A‘e(27,27),B2fl(39) 

C 

DITENSION  INDEX  (352) 

CAa  OPENMSd,  INDEX, 352,0) 

Cifiiiit 

c 

C  READ  INPUT  DATA 

C 

CAa  SETUP 
C 

C  ZERO  THE  FIELDS  INITIAaY. 

C  PROGRAM  USES  A  MASTER  STORAGE  FILE  TO  STORE  FIELDS  AND  CONDUCTIVITY 
C  OF  EACH  LATTICE  POINT.  INDEX  K  GOES  FROM  I  TO  27. 


c 

EX 

K 

c 

EY 

K427 

c 

EZ 

K454(27i2) 

c 

HX 

K461(27f3) 

c 

HY 

K4l0e(27i4) 

c 

HZ 

K4l35(27i5) 

c 

SIS 

K4lG2(27i6) 

C  HXTN2  R4109(2747) 

C  HYTM2  K42l6(27fB) 

C  HZTM2  K4243(27i9) 

C  Cl  K4270(27tl0) 

C  C2  R4297 (27111) 

C  C3  K4324(27fl2) 


ClIHtit 

CiiiiiliftHiDELETE  FROM  HERE  TO  STOP  DELETEi 
C  TO  REMOVE  THE  ZEROING  OF  TAPEI 

C  FOR  CONTINUED  RUNS  BEYOND  TIE 

C  INITIAL  TMAX  ST0PIN6  TIIE 
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CfmifHtltfltHfffiltltltftfifHItfftttltKIlHIH 

DO  ie«  I«i,27 
DO  iM  >1,27 
lee  Ai(i,ji=o.0 

ClIHIIt 

DO  200  K=l,351 

Cfifffff 

200  CALL  URITNSd, Rid, l),7Bfl,K,0) 

C 

CfitfifiHHiSTOP  DELETE  HERE  60  TO  THM= 

C 

C 

C  PRESET  TIE  CONDUCTIVITY  ARRAY 

C 

CALL  RADDC 
C 

C  TIME  STEP  LOOP 

C 

DO  1000  K=1,NMAX 
CALL  APERT 
C 

C  COMPUTE  TIME  CONSTANTS 

C 

CtfiiiiiHiiiFOR  CONTINUED  RUNS  BEYOND  TNAXMMMHit 
C  ADD  THE  STOPPING  TIME  FROM  PREVIOUS 

C  RUNS  TO  THN  AND  TEN  PLUS  1/2  OF  DELTA  T 

C 

C  EXAMPLE 

C 

C  AFTER  FIRST  RUN 
C 

C  THN=DELTi(FL0AT(N)-1.0)d.001B3E-7 

C  TEN=DELTi(FLORT(N)-.5)»I.0O183E-7 

CffHiiiiiiiitftfiiiiiiiiHiififiiiifiiiififiiiiffm 

TIfI=DELTi(a0RT(N)-1.0) 

TEH=DaT»(a00T(N)-.5) 

C 

C  ADVANCE  TIE  H-FIELDS  BY  Z-PLANES 

C 

CALL  HADV 
C 

C  ADVANCE  HE  E-f  lELDS  BY  Z-PLANES 

C 

CALL  ERDV 
1000  CONTINUE 
C 

C  END  OF  TIME  LOOP 

C 

CALL  aOSNSd) 

STOP 

END 

SUBROUTINE  SETUP 
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C  CODE  USES  n  UNIFOm  GRID. 

C 

COmON  mu  NmX,EPST,SIGO,C,TNUO,DELT 

COmON  /SETS/  ICflN,ICMI,ICPl,JCAN,JCMl,JCPl,KCm,KCMI,KCPl 

CONfiN  /SET3/DELX,DELV,DELZ 

Citiitff 

c 

c 

c 

Ckhhi 

C  INPUT  DATA  PROVIDED  BY  USER. 

C 

Cii>»«(i 

C  SET  SIGO^e-B  TO  STOP  ABSORPTION  BC 

DATA  C, EPSO, EPSR, XMUO, SIGQ/3. BEB, 8.  B542E-iB, 1. 8, 1. SSBEE-8, 1. 1/ 

ClHHit 

DATA  ICPl,JCPl,KCPl/27,a7,27/ 

DATA  DELX,DaY,DELZ/.69,.a2,.4S/ 

DATA  THAX/l.E-7/ 

C 

ICAN=ICP1-1 

ICN1=ICAN-1 

JCAN=JCPi-I 

JCN1=JCAN-1 

KCAfH<CPl-i 

KCMI«KCAN-1 

C 

C  VERIFY  TNAT  ICPI,JCP1,KCP1  ARE  ODD. 

C 

IF ( ICAN. EQ. 2i ( ICAN/2) . AND. JCAN. EO. Z* ( JCAN/2) . AND. 

4KCAN.EQ.2«(KCAN/aMG0  TO  1 
PRINT  «,*  THE  NUMBER  OF  GRID  POINTS  MUST  BE  ODD  ’ 

STOP 

1  CONTINUE 

C 

C  BELT  IS  THE  TIME  INCREIENT  (IT  SATISFIES  COURANT  CONDITION). 

C  NMAX  IS  TIE  TOTAL  NUMBER  OF  TIME  INCREMENTS  FIELDS  UOULD  BE  ADVANCED. 

C 

DELT=ANIN1  (DaX,  DELY,  DELZ)  /  (2.  *C) 

NMAX=TMAX/DELT 

EPST=EPSO«EPSR/DELT 

TMUO=DELT/XMUO 

RETURN 

END 

SUBROUTINE  RAOBC 
C 

C  SIG(I,J)  AT  Kl=KtIG2  IS  THE  CONDUCTIVIH  AT  THE  POINT  (X(I),Y(J),Z(K)). 
C  X(I)=(I-.5)iDaX,  Y(JI=(J-.5HDELY,  Z(R)=(K-.5)«DELZ 
C 

CtIHIff 

REAL  RN(27,27),RNP(27,27),THE7A,RC0N 
COMMON  /ARRAYS/  SI6(27,27t,Bl(39) 


n  o  o  c->  r>  rs  r> 


tfli3(27,27),B13(39),AI4(a7,87),BI4(39>,Ci(27,27),B15(39), 
4C2 (27, 27) , B16 (39) , C3 (27, 27) , B17 (39) , niB(27, 27) , BIB (39) , 
4fll9 (27, 27) , B19 (19) , 020 (27, 27) , 220(39) 

(»f1QN  /SET3/  DEL)(,DaV,DELZ 
C 
C 


Cliff ttl 


CONION  /SET!/  NnA7,EPST,SIGQ,C,mJ0,DaT 

COmON  /sn2/  ICflN,IO(l,ICPl,JCAN,JCMi,JCPl,KCnN,KO(l,KCPl 


PRESET  DC  CONDICTIVITY  )ERE  IF  )ODED 


CALCULATE  GEOCTRIC  CONSTANTS  FOR  RA06C 

C  FDDINS  CENTER  OF  SHIFDO  COORDINATE  SYSTEH 
IX0=ICAN/2fl 
1Y0=JCAN/2M 
IZ0=«CAN/2fl 
C 

DO  50  K=1,KCAN 

Kll=Kf270  -V 

K12=«+297 

Ki3=K+324 

c 

CALL  READNS  (1,C1  (1,1), 76B, KID 

CALL  READNS  (1,C2(1,D, 768, K12)  1 

CALL  READHC  (1,C3(1, l),76fl,K13) 

IF  (K.E0.2.0R.K.E0.KCAN)  TIEN  v‘:>; 

C  FfiCE 

KW=IRBS(K-IZ0) 

KNP=KNM 

DO60I=1,I(DN 

DO70J=1,JCAN 

c  :  ■ 

IN=IRBS(I-IX0) 

JN=IABS(J-IY0) 

C 

RN  ( I ,  J)  i^RT  ( (ROAT  ( IN)  iDELX )  f  f  2. 0f  (FLOAT  ( JN)  ilELY)  f  f  2. 0 
♦♦  (FLOAT  (KN)fDELZ)f  12.0) 

RNP  ( I,  J)  :>SORT  ( (FLOAT  ( IN)  f  DELX)  f  f  2. 0f  (FLOAT  ( JN)  f  DRY)  f  f  2. 0 
♦f(aOAT(KNP)fDELZ)ff2.0) 

DETA»(lFIX((RNP(I,J)-RN(I,J))/(DELZ)f0.5))- 

l((RNP(I,J)-RN(I,J))/(DaZ)) 

C 

RCOH=RN(I,J)/(2.0fRNP(I,J)) 

C 

Cl  ( I ,  J) =RCONf  (DETAf  (DCTA-1 . 0) ) 

C2  ( I,  J) =Ra)Kf  2. 0f  ( 1 . 0-  (DETAf  12. 0) )  :■ 

C3(I,J)=RC(]NfDETAf(DETAfl.0) 
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71  CONTINUE 
68  CONTINUE 
C 

ELSE 

C  SIDES  OF  K-PLANES 

KN»inBS(K-IZe) 

C 

icte=icnN-E 

JCN2=JCflN-2 

C  SIDES  OF  1=2  t  ICON 

DO  68  1=2,  ICm,  ICN2 

DO  S8  J=1,JCAN 

C 

iN=inBsn-ix0) 

INP=INM 

JN=IADS(J-IYe) 

C 

RN 1 1,  J)  =SQRT  ( (aOAT  ( IN)  lOaX)  »2. 8^  (ROAT  ( JN)  (DELY)  ii2. 8 
l4(FL0AT(KN)iDELZI*t2.8) 

RNP  ( I,  J)  =SQRT  ( (FLOAT  ( INP)  lOaX)  tf  2. 8^  (FLOAT  (JN)  iDaY)  ttg.  8 
tt(aOAT(KN)iDELZItf2.8l 
C 

T}ETA=(  IFI X  ( (RNP(  I,  J) -RN(  I,  J)  I  /  (DELX)48.  5)  I- 
8((RNP(I,J)-RN(I,J))/(DELX)I 
RC0H=RN(I,J)/(2.8iRNP(I,J)) 

C 

Cl  ( I ,  J) =Ra)Nf  (TltTAi  (TICTA-l .  8) ) 

C2  ( I,  J)  =RCON>2. 8i  ( 1 , 8-  (T)ETAh2.  8) ) 

C3(I,  J)  =RCONi  (TICTA*  (TIETAfl.  8)  I 
C 

98  CONTINUE 
88  CONTINUE 

C  SIDES  AT  J=2  I  JCAN 

DO  180  J=2,JCAN,JCN2 
DO  118  I=1,ICAN 
C 

IN=IRBS(I-IX8) 

JN=IABS(J-1Y8) 

JNP=JN4l 

C 

RN  ( I ,  J)  =SQRT(  (ROOT  (IN)  iDELX )  ii2. 84  (FLOAT  (JN)  iDaY)  ii2. 8 
♦4(aOAT(KN)fDELZ)»»2.8) 

RNP(  I,  J)=SQRT(  (aOAK  IN)  iDELX)  f  12. 84(aOAT(JNP)iDELY)  ii2. 8 
|4(aOAT(KN)iDELZ)fi2.8) 

C 

T)£TA=(IFIX((RNP(I,J)-RN(I,J))/(DaYl48.5))- 

«((RNP(I,J)-RN(I,J))/(DaY)) 

RC(]l1=RN(I,J)/(2.8iRNP(I,J)) 

C 

Cl  ( I ,  J)  =Ra)nt  (T)ETA*  (T)ETA-1 . 8) ) 

C2 ( 1 , J) =RC0Mi2. 8» ( 1. 8- (THETAI42. 8) ) 

C3(  I,  J)  =RCONi  (HETAi  (T)€TA4l.  8) ) 
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n  n  n  n  n  n  n  n  n  onrsojgr* 


111  CONTINUE 
IM  CONTINUE 
C 

END  IF 
C 

CALL  URITMS  (l,Cl(l,n,76a,Kll,l) 
CALLURITHS  (l,C2(l,n,76S,KI8,i) 
CALLURITNS  (1,C3(1, 11,768, K13, 1) 

CONTINUE 


ItfHt 

RETURN 

END 

SUBROUTINE  EAOV 

THIS  SUBROUTINE  ADVANCES  THE  E-FiaOS  TUO  Z-PLANES  AT  A  TINE. 

RECORDS  OF  H-FIELDS  FOR  TUO  SUCCESIVE  Z-PLANES  MUST  BE  KEPT  IN 
ORDER  TO  TAKE  DERIVATIVES  IN  THE  Z-DIRECTION. 

RECORDS  OF  SIS  FOR  TUO  SUCCESIVE  Z-PLA1€S  MUST  BE  KEPT  IN  ORDER  TO  TAKE 
AVERAGES  IN  THE  Z-DIRECTION. 

COMMON  /SETl/fl1AX,EPST,SI60,C,TMU0,DELT 

COrtlON  /SETS/  ICAN,ICM1,ICPI,JCAN,JCM1,JCP1,KCAN,KCM1,KCP1 

COMMON  /SET3/  Dax,Dav,Daz 

COmON  /TIME/  TEN,  TUN 

COmON  /ARRAYS/  EX(a7,27),Bl(39),Ey(27,E7),B2(39),EZ(87,R7), 
4B3(39),HXA(27,S7),B4(39),HXB(27,87),B5<39),HYA(27,S7),B6(39), 
4HYB(27,27),B7(39),HZA(27,87l,Da(39l,HZB(27,27),B9(39), 

4SIGA(27, 871 , Die (39) , SIGD  <87, 87) , B1 1 (39) 

ADVANCE  E-FIELDS  BY  Z-PLANES 
READ  IN  FIRST  H-PLANES 

CALL  REA0MS(I,HXA(1,1),768,68) 

CALL  RERDMS(I,HYA(1,I),768,109) 

DO  FOR  ALL  Z-PLANES 

MNa=B 

)W9=9 

MNie^lB 

DO  IBM  KL00P^1,KCAN,8 

K=«LOOP 

KPl=K4l 

K8=K487 

R3=H454 

K4=K4e8 

K5=«4ie9 

K6H(4135 


n  n  o 


K7>K*162 


C  READ  IN  SECOND  N-PUVES  AND  OLD  E-f  lELDS 
C 

CALL  REAI)NS(i,EX(l,l),7Ca,K) 

CALL  READNS(l,EY(l,n,7Sa,K2) 

CALL  READNS(I,EZ(I,I),768,K3) 

CALL  REA0HS(1,HXB(1,1),7S8,K4) 

CALL  REA0NS(i,HYB(l,l),7Sa,KS) 

CALL  READNS(l,HZA(l,n,7Sa,K6) 

CALL  REnDMS(l,SI6A(l,l),7SB,K7) 

C 

DO  see  >1,JCAN 
JPl=Jtl 

DO  m  I^,  ICAN 

c 

C  SIG(IXD,Y,Zh.5>(SI6(X-DELX,Y,Z)4SIB(X,Y,Z)) 

C  SI6(X,IY8,Z)=.5*(SIG(X,Y-DELY,Z)4SI6(X,Y,Z)I 

C  SIG(X,Y,IZef-.5i(SI6(X,Y,Z-DELZ)4SI6(X,¥,Z)) 

C 

SI62=. 25* (S IGA ( 1, J) 4SIGA ( I-l, J) ) 

A=€PST-SIG2 

B=1,/(EPST*SI62) 

EX(I,J)-(A(EX(I,J)*(HZA(I,JPl)-HZA(I,jn/DELY 
*IHVB(I,J)-HYA(I,J))/DELZ)*B 
m  CONTINUE 

DO  308  I>1,ICAN 
IP1=I*1 

DO  380  J’iGiJCAN 

8162=.  25*  (SIGAI  I,  J-U*SIGO(I,  J) ) 

R=£PST-SIGa 

B=i./IEPST4SI62) 

EYn,J»=(A*EY(I,J)*(HXB(I,J)-HXA(I,J))/DaZ  - 
*(HZA(IP1,J)-4IZA(I,J))/DELXI*B 
380  CONTINUE 

URITE  NEU  EX,EY  PLANES  TO  HASS  STORAK 

CALL  URITHS(I,EX(1,1),768,K,I) 

CALL  URITMS(I,EY(l,l),78a,K2,l) 

IF  (K.EQ.IIGO  TO  450 
DO  400  1=1,  ICAN 
IP1=I*1 

DO  400  J=1,JCAN 
JP1=J*1 

SI62=.25*(SI6AII,J)4SIGBn,JI) 

A=€PST-SIG2 
B=I./(EPST4SIG2) 

EZII,J>=(A*EZ(I,J|4(HYA«IPl,J>-HYAII,J»»/DaX 
4lHXRn,JPl)-NXR(I,JM/DELY)«B 
480  CONTINUE 


CJ  cd  o 


C  URITE  NEU  EZ  RiVC  TO  MASS  STORAGE 
C 

Cflll  URIT)6(1,EZ(1,1>,768,K3,1) 

SAMPLE  POINTS  OF  INTEREST  (ODD  K  PLttCS) 

111  FORMAT  (M,2(1X,1P1E1I.M) 

112  FORMAT  (I4,2(1X,1P1E11.4I,2(/)) 

C 

C 

450  CONTINUE 
K=RP1 
KP1=K*1 
K2=K2tl 
RM3+1 
K4=1<4*1 
K5=«5M 
KS=K6tl 
K7=«7+l 
C 

C  READ  IN  SECOND  OLD  E-PLATES  AND  (EXT  H-PLATES  WER  A  PLANE 
C 

CALL  REAI)M5(l,EX(l,l),7Sa,K) 

CALL  READNS(l,EY(l,l),7Sa,K2) 

CALL  REA0MS(1,EZ(1,1),7SS,K3I 
CALL  READT6(1,HXAI1,1(,768,K4) 

CALL  READNS(1,HYA(I,I),7SB,K5( 

CALL  READMS(l,HZDIl,l),7Sa,KS) 

CALL  READM5(1,SI6B(1,1),768,K7) 

DO  690  M,JCAN 
JP1=J+1 

DO  690  I^.ICAN 

SI62>.2S«(SIGB(I,J)4SIEB(I-I,J)) 

A=EPST-SIG2 

B=1./(EPST4S(62( 

EX(I,JI»(AfEXII,J)4(HZB(I,JPl)-HZB(I,Jn/iaY  - 
4(HYA(I,J(-HYB(I,J))/DELZ)«B 
690  CONTINUE 

DO  799  I«1,ICAN 
IPI-I41 

DO  700  J=2,JCAN 

SI62».2St(SIGB(I,M)4SI6B(I,J)) 

A=€PST-S1G2 

B=1./(EPST4CI62) 

EY(I,J)=(AiEY(I,J)4(HXA(I,J)-HIB(I,J))/DaZ  - 
4(HZBIIPl,J)-HZB(I,JM/DaXliB 
790  CONTINUE 

DO  BOO  U1,ICAN 
IP1>I4| 

DO  BOO  M,JCAN 
JP1»J*1 

SIG2«.2S«ISI6B(I,J)4SIGA(I.J)) 
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A=€PST-6I62 

B»i./(EPST4SI62) 

EZ(I,JN(fliEZ(I,J)4(HYB(IPl,J)-HYB(I,Jn/DaX  - 
4(HXB(I,JPn-HXB(I,J))/DELY)iB 
80e  CONTINUE 
C 

C  URITE  lEU  EX,EY,EZ  PUVCS  TO  HASS  STORAGE 
C 

CALI  URITHS(1,EX(I,1),7SB,K,1) 

CALL  URITNS(l,EY(l,n,>ES,KE,l) 

CALL  HRITHS(1,EZ(1,1),768,K3,1) 

C 

C  SAMPLE  POINTS  OF  INTEREST  (EVEN  K  Pli»£S) 

C 

IF  (K.ED.ie)  URITE  (6,111)  MNa,TEN,EY(17,ll) 

IF  (K.EQ.1B)  URITE  (6,112)  MNte,TEN,EY(17,ll) 

C 

IF  (K.E0.14)  T)EN 

URITE  (6,111)  HH9, TEN, El (15, 13) 
aSE 
END  IF 
C 

1000  CONTINUE 
C 

C  ZERO  TNE  aaTRIC  FiaO  UITHIN  THE  )€TAL  (XMPONENTS  OF  AIRCRAFT. 

C 

CALL  AIRPLN 

RETURN 

END 

SUBROUTINE  HWV 
C 

C  THIS  SUBROUTDE  ADVAfCES  THE  H-FiaOS  TUO  Z-PLANES  AT  A  TIME. 

C  RECORDS  OF  E-FiaOS  FOR  TUO  SUCCESIVE  PLANES  MUST  BE  KEPT  IN  ORDER  TO  TAKE 
C  DERIVATIVES  IN  TIE  Z-DIRECTION. 

CIHIIH 

C  ALSO  MOST  MODIFICATIONS  FOR  T)C  RADIATION 
C  BOUNDARY  CONDITIONS  HERE  PLACED  HERE 

Cffttiif 

c 

COmON  /SETl/  NnAX,EPST,SI60,C,TMUD,DaT 

CO)f«N  /SET2/  ICAN,ICM1,ICP1,JCAN,JCMI,JCP1,KCAN,KCMI,KCP1 

coHtoN  /SET3/  Dax,DaY,Daz 

COmON  /TIME/  TEN,THN 

COHION  /ARRAYS/  EXA(27,27),D1(3S),EXD(27,27),B2(39), 
4EYR(27,27),B3(39),EYB(27, 27), D4(3S),EZA(27,27), 65(39), 
♦EIB(27,27),B6(39),HX (27,27), B7139),HY 127, 27),B8(39), 

4HZ(27,27),B9(39), 

ClIHIII 

4A1B(27,27),BI0(39),A1I(27,27),B11(39),HXTN2(27,27), 612(39), 
4HYTM2(27,27),B13(39),HZTM2(27,27),BU(39),C1(27,27),B15(39), 
4C2(27,27),B16(39),C3(27,27),B17(39),HXKN1  (27,27), BIB(39), 

4HYKM1 (27, 27) , B1 9( 39) ,  HZKMl (27, 27) ,  B20 (39) 


O  o 


REAL  HXTMi (a7,27),HYTHi (27,27),HZTN1  (87,87) 
C 

CtIHIII 


DO  FOR  ALL  Z-PLDNES 
lf(l=l 
(t(2=8 
MN3=3 
»(4=4 
(1(5=5 
Hr(B=6 
(«7=7 
C 

DO  9M  KLD0P=1,KCAN,8 

K=«LOOP 

Krii=«-1 

K8=«*87 

K3=K+54 

K4=Ktai 

K5=KM0a 

K6=Ktl35 


Kia=«+843 

Kll=K^87e 

K18=K+897 

K13=K^3a4 

CALL  READMS  <l,HZrM8(l,l),7a8,Kiai 
CAa  READTIS  (1,C1(1,I), 788, Kill 
CALL  READtC  (1,C8(1,1),768,K18) 

CALL  READMS  (1,C3(1,1),788,K13) 

ClIlfHf 

CALL  REA0N5(1,EXB(1,1),768,K) 

CALL  READMS ( 1, EYB (1,1), 768, K8I 
CALL  READMS 1 1, HZ (1,1 1, 768, K6) 

IF  (K.Ea.l)Ga  TO  358 
Ct»MH» 

K8=KM89 

K9=«*8I6 

CAa  READMS(l,HXTM8(l,l),76e,K8) 

CAa  READMS  (1,HYTM8(1,1),768,K9I 

ClIlHII 

CAa  READM5(1,EZB(1,1),768,K3) 

CAa  READMS(l,HX(t,l),768,K4l 
CALL  READMS(l,HY(l,l),76e,K5) 

DO  888  J=8,JCm 
DO  880  1=1,ICAN 
ClIfllK 

HXTM1(1,J)=HX(1,J) 

CtIHIII 

HX(l,JI=HX(l,J)+T((UOi((EYB(l,J)-EYR(l,J))/DaZ  - 
♦  (EZB(l,J)-EZB(l,J-l))/DaYI 
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£00  CONTINUE 

DO  300  I»i,iail 
1P1=I+1 

DO  300  J=i,JCflN 


HYTH1UP1,J)=HY(1P1,J) 

Ciiifft# 

HY(IPl,J)=«Y(lPI,J)*TNUO*((EZB(IPl,JI-E2Bn,J»»/DELX  - 
+  (EXB(IPl,JI-EXfl(IPl,J))/DaZ) 

300  CONTINUE 
C 

c  flPPRoxinnTES  h  fields  in  plane  perpendicular  to  wire 

C  BY  H=I/eiPIiR. 

C  SOURCE  FUNCTION, ODD  K  REFERENCE  PLAf£S,  HY 

C  CURRENT  PULSE  HOVING  IN  X  DIRECTION  ATTADES  TO  TTE  NOSE  OF  TIE  AIRPLANE 
C  AT  (X,Y,Z)>(3,9,!5),  EXITS  BY  THE  TAIL  AT  (X,Y,Z>>=(£5, 10, 15). 

C  AT  Tl)t  T=0.0  PULSE  IS  LOCATED  AT  X0t3). 

C 

1F(K.NE.15)G0  TO  1 
DO  10  1=2, 3,1 

TPULSE=THW-. 69i (FLOAT ( I ) -3. ) /C 
1F(TPULSE.LE.0.0)GO  TO  10 
C 

X=EXP(-3.5E0MTPULSE)-EXP(-2.625E7iTPULCE) 

C 

HYS=7.07E3iX 
HY(1,9I=-HY9 
10  CONTINUE 

DO  30  1=25,26,1 

TPULSE=TIN-.  691  (ROAT  ( 1 1  -3. 0)  /C 
IF(TPULSE,LE.0.0)GO  TO  30 
l=EXP(-3. 5E4iTPULSEI-EXP(-2. 625E7fTPULBE) 

HYS=7.07E3»X 
HY(I,10)=-HYS 
30  CONTINUE 

1  CONTINUE 

CfllHIf 

C  SIDES  WIEN  1=1  tlCPl 

ICM2=ICAN-2 
DO  310  1=2, ICON, ICH2 
DO  320  J=1,JCAN 
IF  (I.E0.2)UEN 
IB=1 
ELSE 

IB=ICANM 

EM)  IF 

HY  ( IB,  J)  =C1  ( I,  J)  iHYTTB  ( I,  J) +C2  ( I,  J)  tHYTNl  ( I,  J) +031 1,  J)  iHY  ( I ,  J) 
HYTN2(I,J)=HYTN1(I,J) 

IF  (J.BE.2)  THEN 

HZ(IB,J)=Cl(I,J)«HZTM2(I,J)tC2(I,J)«HZTHl(I,J)4C3(I,J)iHZ(I,J) 

HITN2(I,J)=HZTH11I,J) 


END  IF 

3SD  CONTINUE 
310  CONTINUE 
C  SIDES  UlCN  >1  I  JCPl 
JCM2=JCflN-E 
DO  330  J=2,JCnN,JC)12 
DO  340  I=!,ICAN 
IF  (J.EO.S)  THEN 
JB=I 
ELSE 

JB=JCflNM 

END  IF 

HX(I,JB)=Cl(I,J)»HXTK2(I,JUC2(I,JI»HXTW(I,J)+C3n,J)»Hni,J) 

HXTI12(I,J)=HXTMi(I,J) 

IF  (I.6E.8)  TIEN 

HZ(I,JB)=Cl(I,J)fHZTN2(I,J)+C2(I,JliHZTNl(I,J)4C3(I,J)iHZ(I,J) 

HZTK2(I,J)=HZTM1(I,J> 

ELSE 
END  IF 

340  CONTINUE 
330  CONTINUE 

HZ(I,1I=HZ(2,1) 

HZ(l,a7)=HZ<2,27» 

HZ(27,l»=HZ(27,2l 

HZ(27,*7»=HZ(27,26) 

Ciitiii* 

c 

%  C  WRITE  NEW  HX,W  PlflMES  TO  NftSS  STORAGE 

C 

Cfla  HRITNS(I,HX(l,l),760,K4,n 
CALL  HRITNSU,HY(l,n,768,K5,l» 

Ciitiiit 

CALL  HRITNS(l,HZ(I,l),76fl,K6,l) 

CALL  HRIIMS(I,HXTK2(1,1),760,KS,I) 

CALL  HRITNS(l,HYTM2(l,U,760,K9,l) 

CAU  WRITNS(l,HZTN2ll,l),7M,K10,I) 

CfllHH 
350  CONTINUE 

DO  400  I=I,ICN1 
IP1=I+1 

DO  400  J=1,JCN1 
JPI=Jtl 

CimiH 

HZTN1(IPI,JP1)=HZ(IP1,JP1) 

ClIHIII 

HZ(IPl,JPI)=HZ(IPl,JPl)tTW)0«((EXB(IPl,JPl)-EXB(IPl,J)»/DaY  - 
♦  (EYB(IP1,JP1)-EYB(I,JP1))/Dax» 

400  CONTINUE 

Ctftifit 

C  SIDES  IM  I  ICPt 
iai2=ICAN-2 
DO  610  1=2, ICAN, ICN2 
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DO  £20  J=2,JCAN 
IF  (I.ED.E)  T>€N 
IB=i 
ELSE 

IB^lCPNfl 

END  IF 
C 

HZ(IB,J)<^1(I,J)iHZTN2(I,J)H:2(I,J)<HZTN1II,J)4C3(I,J)iHZ(I,J) 

HZTK2(I,J)=HZTHI(I,J) 

£20  CONTINUE 
£10  CONTINUE 
JCM2=JCflN-2 

C  SIDES  J=1  I  JCPi 
DO  £30  J=2,JCftN,JCN2 
DO  £A0  1=2, ICON 
IF  (J.E0.2)  THEN 
JB=1 
asE 

JB=JCflNM 

END  IF 

HZ(I,JB)=CI(I,JHHZTN2(I,J)+C2(l,JHHZTNin,J)4C3n,J)tHZ(I,J) 

HJTM2(I,J)=HZTN1(I,J) 

£40  CONTINUE 
£30  CONTINUE 

HZ(1,1)=HZ(2,1) 

HZ(I,27)=HZ(2,27) 

HZ  127, 1 1 =HZ (27, 2) 

HI(27,27I=HI(27,2£» 

C  WRITE  NEW  HZ  PLATC  TO  MOSS  STORAGE 
C 

COLL  URITMS(l,HZ(l,l),7£e,K£,l) 

COLL  URITMS(l,HZTH2ll,l),7£e,K10,l) 

C»fnm 

C 

C  SAMPLE  POINTS  OF  INTEREST  (ODD  K  PLANES) 

C 

IF  (K.EA.13)  URITE(£,111)  )V(2,THN,HZ(4,10) 

IF  (K.EQ.13)  HRITE(£,II1)  MN3,THN,HZ(24, II) 

IF  (K.EQ.23)  URITE(£,111)  MN7,THN,H)((18, 10) 

C 

111  FORMAT  (I4,2(1X,1P1E11.4)) 

C 

K=K+I 

KNI=K-1 

K2=K2M 

K3=«3n 

K4=K4+1 

K5=H5+1 

K£=H£+I 

Cii«i«ii 

K8=«Ma9 

K9=)(>21£ 


108 


• 


Kie=K«243 

Kll=Kt27g 

K12=R+E97 

R13=Kt324 

CALL  REAI»<5(l,HXTM2n,l),7SB,K8) 

CAIl  REnDI1S(I,HYTI12(l,I),7&8,K9) 

CALL  READNS(1,HZTI12(I,1),768,K1  ) 

CALL  READI1SU,CI(1,1),768,K1I) 

CALL  READNS(1,C2(1,I),7C8,K12) 

CALL  READMS(1,C3(1,1),76B,K13) 

C  FACE 

IF  (K.EQ.2.0R.K.Ea.KCAN)  THEN 
IF  (K.E0.2)  THEN 

KMIXrKtee 

KNlY=Kt|07 

RNIZ=«*134 

CALL  READMS  (l,HXKm(l,l),768,KmX) 

CALL  READNS  (1,HYKMI(I,1),768,KM1Y) 

CALL  READNS  ll,HZKMl(l,l),76a,KNlZ) 

ELSE 

KNlX=Kt82 

KHIY=«tl09 

KN1Z=KM36 

CALL  READNS  (l.HXKMUl.n.TEB.KMUI 
CALL  READNS  (I,HYKNl(l,n,7M,KNlY) 

CALL  READNS  (1,HZKN1(1,1),768,KMIZI 
END  IF 
ELSE 
END  IF 

Cffiffii 

c 

C  BEAD  IN  NEXT  PLAINS  OF  E  AND  H 
C 

CALL  READNS(l,EXA(l,l),768,K) 

CALL  READNSIl,EYA(l,l),768,K2) 

CALL  READNS(1,EZA(I,1),7S8,K3) 

CALL  REftDNS(l,HX(l,U,768,K4» 

CALL  RERDNS(I,HY(I,I),768,K5) 

CALL  READNS(1,HZ(1,U,7S8,K6) 

DO  680  J=1,JCN1 
JP1=J+1 

DO  600  I=1,ICAN 

ClKHtf 

HXTNl(I,JPn:4IX(I,JPl) 

Cdiiiif 

HX(I,JPl)=«X(I,JPl)*TNUO»((EYA(I,JPn-EYB(I,JPn)/DELZ 
♦  (EZA(I,JPU-EZA(I,JM/DELYI 
600  CONTINUE 

DO  700  H1,ICN1 
IPM+1 

DO  700  J=1,JCW 

C«i»ti«« 
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CilHttf 

HY(IPl,J)^1Y(IPt,J)4TflU0f((EZn(IPl,J)-EZfl(I,J))/DaX  - 
«(EXn(IPl,J)-€XB(IPl,J))/DELZt 
700  CONTINUE 
C 

C  SOURCE  FUNCTION,  EVEN  K  REFERENCE  PLiVES,  KY 
C 

IF  (K.NE.  14)00  TO  2 
DO  20  1=2, 3,1 

TPULSE=THN-. 69i (FLOAT ( I ) -3. ) /C 
IF  (TPULSE.LE.  0.0)00  TO  20 

X=EXP(-3. 5E4ITPULSE) -EXP (-2. 625E7tTPULSE) 

HYS=7.07E3iX 
HY(l,9)=tHYS 
20  CONTINUE 

00  40  1=25,26,1 

TPULSE=T)(M-.  69i  (aOAT  ( I )  -3. )  /C 
IF  (TPULSE.LE.0.0)6O  TO  40 

X-EXP(-3.5E4iTPULS£)-EXP(-2.625E7#TPULSE) 

HYB=7.07E3tX 
HY(1,10)=HYS 
40  CONTINUE 

2  CONTINUE 

DO  B00  1=1,ICN1 
IPl=l+l 

DO  B00  J=1,JCN1 
JPI=J*1 

ClItKfl 

HZTN1(IP1,JP1)=HZ(IP1,JP1) 

Ciiiitt* 

HZ(IPl,JPl)=HZ(IPl,JPl)+T)1U0i((EXn(IPI,JPI)-EXA(IPl,J))/DaY 

♦(EYR(IPl,JPl)-EYn(I,JPl))/DaX) 

600  CONTINUE 

C 

C 

C  SOURCE  FUNCTION,  EVEN  K  REFERENCE  PLANE,  HZ 

IFIK.NE. 14)00  TO  4 
DO  60  1=2,3, 1 

TPaSE=THN-.  69i  (FLOAT  ( 1 )  -3. )  /C 
IF(TPaSE.LE.0.0)OO  TO  60 

X=  (EXP  (-3.  K4fTPULSE)  -EXP  (-2. 625E7«TPU.SE) ) 

HZS=1.447E4fX 
HZ(1,9)=-HZB 
HZ(I,10)=HZS 
60  CONTINUE 

DO  60  1=25,26,1 

TPULSE=THN-. 69i (FLOAT ( I ) -3. ) /C 
IF(TPULSE.LE.0.0)6O  TO  60 

X= (EXP  (-3. 5E4fTPUL5E)  -EXP  (-2. 625E7«TPl)LSE) ) 

HZS=1.447E4iX 

HZ(I,I0)=-HZS 


HZ(I,11)=HZS 
80  CONTINUE 

4  CONTINUE 


C  FACE 

IF  (K.EO.a.OR.K.EQ.KCAN)TTCN 
DO  210  1=2, ICON 
DO  220  J=2,JCm 

HXKI1I(I,J)sCl(I,J)*HXTM2(I,J|iC2(I,J)«HXTNl(I,J)4C3(I,J)iHZ(I,J) 

HXTN2n,J)=HXTN!(I,J) 

HYKm(I,J)=Cl(I,J)»HYTH2(I,J)+C2(I,J)»HYTNin,J)4C3n,JUHY(l,J) 

HYTN2{I,J)=HYTM1(I,JI 

HZKNHI,J)^l(I,J)»HZTH2(I,J)4C2(I,J)iHZTHiU,J)4C3(l,J><HZ(I,J) 

HZTM2II,J)=HZTN!(I,J) 

220  CONTINUE 
210  CONTINUE 
C  EDGES  OF  FACES 
lCN2=ICAN-2 
DO  410  1=2, ICAN, ICN2 
DO  420  J=2,JCAN 
IF  (I.EQ.2)  THEN 
IB=1 

ELSE 

IB=ICAN4l 

END  IF 

HXKN1(IB,J)=NXKM1(I,J) 

420  CONTINUE 
410  CONTINUE 
JCM2=JCAN-2 
DO  430  J=2,JCAN,JC«2 
DO  440  1=2, ICAN 
IF  (J.EQ.2)  THEN 
JB=1 
ELSE 

JB=JCAN4l 

END  IF 

HYKM1(I,JB)=HYKH1(I,J) 

440  CONTINUE 
430  CONTINUE 

cna  URITMS(l,HXKHl(l,l),768,KMlX,n 
CAa  URITN5(l,HYKNin,n,768,KMlY,ll 
CAa  URITNS(l,HZKNl(l,l),76a,KHlZ,n 
asE 
END  IF 
IC)12=ICAN-2 

C  SIDED  1=  1  I  ICPl 
DO  230  1=2,  ICAN,  ICI12 
DO  240  J=1,JCAN 
IF  (I.E0.2)  TTEN 
IB=1 
asE 

IB=ICAN4l 
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END  IF 

HY(IB,J)<l(I,J)«HYT112(I,J)4Ca(I,J)iH)mil(I,J)4C3(I,J)iHY(I,JI 

KYTK2(I,J)<4{YTNl(I,J) 

IF  U.GE.2)  TIEN 

HZ(IB,JXI(I,J)iHZTM2(I,J)H2(I,J)tHZTNl(l,J)4C3(l,J)fHZ(I,J) 

HZTH2(I,J)mZTHI(I,J) 

ELSE 
END  IF 

2A0  CONTINUE 
230  CONTINUE 
JCH2=JCflN-2 
DO  250  J=2,JCnN,JCN2 
DO  260  i»i,icm 
IF  (J.EQ.2)  TIEN 
JB=1 

ELSE 

JB«JCflNM 

END  IF 

H)((I,JB)>CI(I,J)«HXTN2(I,J)H:2(I,J)iHITNi(I,J)4C3(I,J)«Hni,J> 

HXTN2II,J)=HXTHt(I,J) 

IF  (I.GE.2)  TIEN 

HZn,JBKl(I,J)<HZTN2(I,J)H:2(I,J)*HZTNl(I,J)4C3(I,J)fHZ(I,J) 

HZTN2(I,J)=NZTN1(I,J) 

ELSE 
END  IF 

260  CONTINUE 
250  CONTINUE 

HZ(l,n=NZ(2,I) 

HZ (I, 27) =HZ 12, 271 
HZ(27,1)=HZ(27,2) 

HZ (27, 271 =HZ (27, 26) 

COLLURITHS  (l,HXTN2(l,l),76fl,KS,  1) 

CflLLNRITHS  (I,HYTH2(1, 1),76S,K9, 1) 

COLLHRITNS  (1,HZTN2(1,1),768,K10,I) 

CfIKIff 

CALL  URITHS(1,HX(I,1),76B,KS1) 

COLL  URITNS(I,HY(I,1),76S,KS,1) 

COLL  URITNS(1,HZ(1,1),76B,K6,1) 

C 

C  SONPLE  POINTS  OF  INTEREST  (EVEN  K  PIAIES) 

C 

IF  (K.EO.U)  THEN 

URITE  (6,111)  HN4,TIII,HZ(I0,4) 

WRITE  (6,111)  NN5,T)M,HZ(I9,4) 

URIIE  (6,111)  HN6,TIt(,HX(22,lS) 

ELSE 
END  IF 
C 

IF  (K.E0.6)  URITE(6,111)  If(l,THN,HX(ie,lD) 

C 

900  CONTINUE 
RETURN 


END 

SUBROUTirC  APERT 

RETURN 

END 

SUBROUTINE  AIRPLN 
C 

C  ZEROES  THE  aECTRIC  FiaDS  WITHIN  THE  Fit  AIRCRAFT  AND 
C  TANGENTIAL  aECTRIC  FiaO  CONPONENTS  AT  TIE  SURFACE 
C 

COmON  /ARRAYS/  FLD  (S7,a7),Bl  (39) 

C 

C  EX  ZEROED 

C 

DO  lem  K=3,2S 

CALL  READNS  ll,FLD(l,l),768,K) 

C  WING 

IF  (K.NE.3)  GO  TO  40 


40 

FLD  (10, 9)  =0.0 

FLD  (19, 9)  =0.0 

IF  IK.NE.4)  GO  TO  50 

41 

DO  41  1=17,19 

FLD(  1, 91=0.0 

50 

IF  (K.LT.S.OR.K.GT.G)  GO  TO  70 

51 

DO  51  1=16,19 

FLO(I,9)»0.0 

70 

IF  (K.LT.7.0R.K.GT.a)  GO  TO  90 

71 

DO  71  1=15,19 

DO  71  J=9,10 
aO(I,J)=0.0 

90 

ao  (23, 9)  =0.0 
aD(24, 91=0.0 

IF  (K.NE.9)  GO  TO  100 

91 

DO  91  1=14,19 

DO  91  J=9, 10 

FLD(I,J)=0.0 

92 

DO  92  1=22,24 
aO(l,9)=0.0 

100 

IF  IK.NE.10)  GO  TO  110 

101 

DO  101  1=11,19 

DO  101  J=9, 10 
aD(I,J)=0.0 

102 

DO  102  1=21,24 
aD(I,9)=0.0 

110 

IF  (K.NE.ll)  GO  TO  120 

111 

DO  111  1=11,23 

DO  111  J=9,10 
aD(I,J)=0.0 

112 

DO  112  1=11,22 
aD(I,B)=0.0 

C 

WAIN  FUSaAGE 

120 

IF  (K.LT.12.0R.K.GT.1S)  GO  TO  126 

DO  121  1=10,19 

113 


u  u  u 


DO  lat  J=5, 18 

181  FLD(I,J)=I.I 

FORUARD  FUSELAGE 

DO  188  1=6,9 
DO  182  J=B, le 

182  FLDn,J)=e.l 
DO  183  1=6,9 
DO  183  J=ll,18 

183  aD(I,J)=6.6 
FU)(7,ll)=0.t 

C  REAR  FUSELAGE 
DO  124  1=80,22 

184  FLD(1,6)=0.0 
DO  185  1=20,83 
DO  185  J=7, 12 

185  FLD(I,J»=0.0 

186  IF  (K.NE.12.0R.K.1E.15)  GO  TO  130 
DO  127  1=19,80 

DO  127  J=3,4 
127  aO  1I,J)=0.0 
aD(20,5)=0.0 

130  IF  (K.LT.13.0R.K.GT.14)  GO  TO  135 
aD(24, 91=0.0 

aD(24, 101=0.0 
C  NOSE  AND  COCKPIT 
DO  131  1=4,5 
DO  131  J=8,9 

131  aO(l,Jl=0.0 
aD(5, 101=0.0 
aD(6, 111=0.0 
aD(7, 181=0.0 
DO  132  1=8,11 
DO  132  J=13,14 

138  aon,Ji=0.0 
aD(I2, 131=0.0 
FLO(  13, 131=0.0 
FLD (9, 151=0.0 
C 

C  RIGHT  VERTICAL  STADALIZER 
C 

135  IF  (K.NE. 131  GO  TO  140 
DO  136  1=19,83 

DO  136  3=13,14 

136  aD(I,Jl=0.0 
aD(17, 131=0.0 
a0(18, 131=0.0 

140  IF  (K.NE. 141  GO  TO  160 
DO  141  1=18,83 

DO  141  3=13,14 

141  aOd, 31=0.0 
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FLO(17,23)=0.I 
DO  142  1=20,23 
DO  U2  >15,19 

142  FLD(I,J)=I.I 
FLD(2e,  15)=e.i 
FLO(28,16)=e.O 
DO  143  1=22,23 
DO  143  3=20,22 

143  FLDd, 31=0.0 
FLD(24, 211=8.0 
FID(24,221=0.0 
DO  144  1=23,24 
DO  144  3=23,24 

144  FLD(1,31=0.0 

160  IF  (K.NE. 161  GO  TO  170 
DO  161  1=11,22 

DO  161  3=0,10 

161  FLD(I,31=0.O 
aD(23, 91=0.0 
FLD  (23, 101=0.0 

170  IF  (K.NE. 171  GO  TO  100 

C  UIN6 

DO  171  1=13,19 
DO  171  3=9, 10 

171  FLO(I,31=e.O 
C 

C  RIGHT  HORIZONTAL  STABOLIZER 

C 

DO  172  1=21,24 

172  FLD(I,91=0.0 

100  IF  (K.l€. 101  GO  TO  190 
DO  101  1=14,19 

DO  101  3=9, 10 

101  FLD(I,3)=0.0 
DO  102  1=22,24 

102  aD(I,9)=0.0 

190  IF  (K.LT.19.0R.K.6T.20)  GO 
DO  191  1=15,19 

DO  191  3=9,10 

191  aD(I,3)=0.0 
aD  (23, 9)  =0.0 
aD  (24, 9)  =0.0 

210  IF  (K.LT.21.0R.K.6T.22)  GO 
DO  211  1=16,19 

211  aD(l,9)=0.0 

230  IF  (K.NE. 231  GO  TO  240 
DO  231  1=17,19 

231  FLD(I,9)=0.0 

240  IF  (K.NE.24)  GO  TO  1000 
aOdO,  91=0.0 
aD(19,9)=0.0 

C  STORE  EX  FIELDS  BACK  IN  TO 


lOM  CALL  URlT11S(l,aD(l,U,«6B,K,l) 

C  ZERO  EY 

DO  aecn  k=>3,24 
L=K427 

GALL  REWNS  (1,FID(1,1),76B,L) 

IF  (K.LT.7.0R.K.GT.ie)  60  TO  lUB 
DO  1081  1^4,19 
10B1  FLD(I,1 01=0.0 

IF  (K.NE.9)  GO  TO  1100 
n.D(13, 101=0.0 

1100  IF  (K.NE. 10)  GO  TO  1110 
FLD(12, 10)=0.0 
FLD(  13, 101=0.0 

1110  IF  IK.NE.ll)  GO  TO  1120 

DO  nil  1=10,22 
DO  nil  J=9,10 
nil  FLD(I,J)=0.B 

FLD(23, 10)  =0.0 
C  MAIN  FUSaflGE 

1120  IF  (K.LT.12.0R.K.GT.15)  GO  TO  1130 
DO  1121  1=9,19 

DO  1121  J=G,12 

1121  FLD(I,J)=0.0 

C  FUD  FUSELAGE 

DO  1122  1=S,B 
DO  1122  J=9,10 

1122  aO(I,J)=0.O 
DO  1123  1=7, B 
DO  1123  J=n,12 

1123  FLDII,J)=0.0 

aD(6,n)=O.0 
C  REAR  FUSaAGE 

DO  1124  1=20,22 
DO  1124  J=7,12 

1124  aD(I,J)=0.0 
DO  1125  J=B,12 

1125  aD(23,J)=0.0 

1130  IF  <K.LT.13.0R.K.6T.14)  GO  TO  1126 

FLD(24,1O)=0.O 

1126  IF  iK.)C.I2.0R.K.NE.15)  GOTO  1131 

DO  1127  I=1B,20 

DO  1127  J=4,5 

1127  FLD  (I,J)=0.B 

FID  (20, 6)  =0.0 

1131  IF  (K.LT.I3.0R.K.GT.14)  GO  TO  1135 

aD(3,9)=0.0 

FLD(4,9)=0.0 
aDI4,10)=0.O 
FLD  (5, 11)  =0.0 
aD(6,12)=0.0 
DO  1132  1=7,13 

1132  FU)II,13)=0.O 
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FL0(I,14)«fl.a 

FL0(8,iS)=e.e 

FLD(9,is)=e.e 

IF  IK.NE.13)  GO  TO 
DO  1136  I>16,83 
FIJ)(I,13)=9.8 
DO  1137  I=ia,23 
FLD(I,14)<^.I 
IF  (K.1C.14)  GO  TO  1160 
DO  lUl  1=22,23 
DO  1141  J=13,24 
FU)II,J)=0.l 
DO  1142  J=21,24 
FLD(24,J)=0.O 
DO  1143  J=20,22 
FU)(21,J)=0.O 
DO  1144  J=17,19 
DO  1144  1=20,21 
aD(I,JI=0.0 
DO  1145  1=19,21 
DO  1145  J=15,16 
aD(I,J)=0.8 
DO  1146  1=17,21 
aD(I,14)=0.0 
DO  1147  1=16,21 
aD(I,13)=0.0 
IF  (K.NE.16)  60  TO  1! 70 
DO  1161  1=10,22 
DO  !161  J=9,10 
aDn,J)=0.0 
aD(23, 101-0.0 
UIN6 

IF  (K.NE.17)  GO  TO  1160 
DO  1171  1=12,19 
FLD(I,10)=0.0 
IF  (K.NE.18)  GO  TO  1190 
DO  1161  1=13,19 
FU)n,10)=0.0 

IF  (K.LT.19.OR.K.6T.20)  GO  TO  i 
DO  1191  1=14,19 
FLD(I,10)=0.0 

GTORE  EY  BflCK  TO  MRSS  STORflGE 
CALL  HRITHS(1,FU)(1, 11,768,1,11 
EZ  ZEROED 
DO  3000  K=3,24 
L=Kt54 

CflLL  REfiD11S(l,aD{l,l),768,L) 

UIN6  TIP 

IF  (K.NE.4)  GO  TO  2050 
DO  2041  1=17,19 

aD(I,9)=0.0 
IF  (K.1C.S)  GO  TO  2060 


DO  eesi  i>i6,i9 
S0S1  FLD(I,9):=e.l 

C  PLANE 

EOED  IF  (K.LT.B.OR.K.GT.7i  GO  TO  2000 
DO  2061  1=15, 19 
2061  FLDII, 91=0.0 

C  WING 

20B0  IF  (K.NE.B)  GO  TO  2090 

DO  2081  J=9,10 
DO  2081  I=1A,19 
20(1  FLD(I,J)=0.D 

C 

C  RIGHT  HAND  GTABALIZER 

C 

DO  2082  1=22,24 
2082  RDd, 91=0.0 

2030  IF  (K.NE.91  GO  TO  2100 

DO  2091  1=14,19 
DO  2091  J=9, 10 

2091  FLD(I,J1=O.0 
DO  2092  1=22,24 

2092  aD(l,91=0.0 

2100  IF  (K.IC.  101  GO  TO  2110 

DO  2101  1=13,19 

DO  2101  J=9,10 

2101  FLD(I,J1=0.0 
DO  2102  1=21,24 

2102  FL0<I,91=0.0 

2110  IF  (K.NE. Ill  GO  TO  2120 

DO  2111  1=12,24 

2111  FLD(I,91=0.0 
DO  2112  1=12,19 

2112  FIO(1,101=O.O 

2120  IF  (K.NE. 121  GO  TO  2130 

DO  2121  1=10,22 

DO  2121  J=0,1O 

2121  FLDII,J1=O.0 
FLD(23,91=0.O 
FU}(23,1O1=O.0 

2130  IF  IK.LT.13.0R.K.GT.1S1  GO  TO  2140 

DO  2131  1=9,19 

DO  2131  J=5, 12 

2131  aO(I,Jl=0.0 

DO  2132  1=5,8 

DO  2132  J=B,10 

2132  FLDII, J1=0.0 

DO  2133  1=7,0 

DO  2133  J=ll,12 

2133  aD(I,Jl=0.0 

FLD (6, 111=0.0 
DO  2134  1=20,22 

DO  2134  J=6,12 


V-X-vj: 


213*  aD(I,J)=B.i 

DO  2135  J=7,12 
2135  FLD(23,J)»«.0 

2140  IF  (K.NE.14)  GO  TO  21B0 

FID  (24, 9)  =0.0 

aD(24, 101=0.0 
DO  2141  1=3,4 
DO  2141  J=a,9 

2141  aD(l,J)=0.0 
aD(4, 101=0.0 
aD(5, 111=0.0 
FLD (6, 121=0.0 

DO  2142  1=7,11 
DO  2142  J=13,14 

2142  aD(I,Jl=8.0 
aD(12, 131=0.0 
FU1(13, 131=0.0 
aD(8, 151=0.0 
FIJI  (9, 151=0.0 

DO  2143  1=16,23 

2143  aD(1, 131=0.0 
DO  2144  1=18,23 

2144  aDd,  141=0.0 

2160  IF  IK. NE. 161  GO  TO  2170 
DO  2161  1=10,22 

DO  2161  J=8,10 

2161  aDII,Jl=0.0 
aO  (23, 91  =0.0 
FLD  (23, 101=0.0 

2170  IF  (K.NE. 171  GO  TO  2180 
DO  2171  1=12,24 

2171  FLDIl, 91=0.0 
DO  2172  1=12,19 

2172  FLDd,  101=0.0 

2180  IF  (K.NE. 181  GO  TO  2190 
DO  2181  1=13,19 

DO  2181  J=9,10 

2181  aD(I,Jl=0.0 
DO  2182  1=21,24 

2182  aDd, 91=0.0 

2190  IF  (K.LT.19.OR.K.GT.201  GO  TO  2210 
DO  2191  1=14,19 

DO  2191  J=9,10 

2191  FLDd,Jl=0.0 
DO  2192  1=22,24 

2192  FLDd,91=0.0 

2210  IF  IK.LT.81.0R.K.GT.221  GO  TO  2230 
DO  2211  1=15,19 

2211  FLDd, 91=0.0 

2230  IF(K.  IE. 231  GO  TO  2240 

DO  2231  1=16,19 

2231  FLDd, 91=8.0 
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8240  IF  (K.NE.24)  GO  TO  3009 
DO  2241  I>17, 19 

,  2241  F1J)(I,3)=0.I 

I  c 

C  STORE  EZ  BflCK  TO  MSS  STORflGE 
C 

3000  CALL  WRITHS  n,FLD(l,l),7M,L,l» 

RETURN 

■  END 

I  END  OF  FILE 


it 


1 


I 


i 


I 


I 


I 


120 


Appendix  C 


Sensor  Plots 


This  appendix  contains  graphs  of  the  ten  sensors  for 
both  the  radiation  and  absorption  boundary  conditions.  Each 
page  contains  one  graph  for  the  absorption  boundary 
conditions  and  one  graph  for  the  radiation  boundary 
condition.  Each  graph  is  made  up  of  1000  data  points 
calculated  from  the  corresponding  programs  in  Appendix  A  or 
B.  The  points  were  connected  with  a  natural  cubic  spline. 
The  only  differences  in  the  two  codes  are  the  boundary 
conditions  as  annotated  in  the  program.  Both  plots  on  the 
same  page  are  for  identical  time  periods  so  that  side-by- 
side  comparisons  could  be  made.  The  programs  were  run  for 
results  out  to  2.5  microseconds.  The  graphs  in  this 
appendix  only  reflect  the  first  1.0  microseconds.  The 
graphs  beyond  1.0  microseconds  were  only  continuations  of 
the  same  trends  as  previous  graphs,  so  they  were  deleted  for 
brevity.  Also,  Absorption  and  Radiation  Boundary 
Conditions  are  abbreviated  ABC  and  RBC  respectively  on  the 


Figure  27.  Sensor  One,  H-field,  Right  Wing,  H  (18,10,6) 


Figure  29.  Sensor  One,  H-field,  Right  Wing,  H  (18,10,6) 
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Figure  30.  Sensor  One,  H-field,  Right  Wing,  H  (18,10,6) 


Figure  32.  Sensor  Two,  H-field,  Nose 


Figure  35.  Sensor  Three,  H-field,  Engine  Burner  Ca: 
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Sensor  Three,  H-field,  Engine  Burner  Ca: 


Three,  H-field,  Engine  Burner  Ca: 


Sensor  Four,  H-field,  Forward  Fuselage  Bottom 


-asc 


Figure  42.  Sensor  Four,  H-field,  Forward  Fuselage  Bottom 


ao.oon 


Five,  H-field,  Rear  Fuselage  Bottom 


Sensor  Five,  H-field,  Rear  Fuselage  Bottom 


Sensor  Five,  H-field,  Rear  Fuselage  Bottom 


Figure  46.  Sensor  Five,  H-field,  Rear  Fuselage  Bottom,  H  (19,4,14) 


Sensor  Six,  H-field,  Vertical  Stabilize 


Figure  53*  Sensor  Seven,  H-field,  Left  Wing 


Sensor  Eight,  E-field,  Right  Wing,  E  (17,11,10) 


Eight,  E-field,  Right  Wing, 


Figure  58.  Sensor  Eight,  E-field,  Right  Wing,  E  (17,11,10) 


Figure  59.  Sensor  Nine,  E-field,  Middle  Fuselage  Top,  E  (15,13,14) 


Sensor  Nine,  E-field,  Middle  Fuselage  Top 


Figure  62.  Sensor  Nine,  E-field,  Middle  Fuselage  Top,  E  (15,13,14) 


Sensor  Ten,  E-field,  Left  Wing 


qure  64.  Sensor  Ten,  E-field,  Left  Wing 
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Sensor  Ten,  E-field,  Left  Wing 
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